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v  Abstract 

Computational  modules  early  in  the  human  vision  system  typically  generate  sparse  information 
about  the  shapes  of  visible  surfaces  in  the  scene.  Moreover,  visual  processes  such  as  stereopsis  can 
provide  such  information  at  a  number  of  levels  spanning  a  range  of  resolutions.  In  this  paper,  we 
extemf-this  multi-level  structure^to  encompass  the  subsequent  task  of  reconstructing  full  surface 
descriptions  from  the  sparse  information^The  mathematical  development  proceeds  in  three  step*. 
First,  the  surface  most  consistent  with  the  sparse  constraints  is  characterized  as  the  solution  to  an 
optimal  approximation  problem  which  is  posed  as  a  variational  principle  describing  the  constrained 
equilibrium  state  of  a  thin  flexible  plate.  Second,  local,  finite  element  representations  of  surfaces 
are  introduced  and,  by  applying  the  finite  element  method,  the  continuous  variational  principle 
is  transformed  into  a  discrete  problem  in  the  form  of  a  large  system  of  linear  algebraic  equations 
whose  solution  is  computable  by  local- support,  cooperative  mechanisms.  Third,  to  exploit  the 
information  available  at  each  level  of  resolution,  a  hierarchy  of  discrete  problems  is  formulated 
and  a  highly  efficient  multi-level  algorithm,  involving  both  intra-level  relaxation  processes  and 
bi-directional  inter-level  local  interpolation  processes  is  applied  to  their  simultaneous  solution. 
Examples  of  the  generation  of  hierarchies  of  surface  representations  from  stereo  constraint*  are 
given.  Finally,  the  basic  surface  approximation  problem  is  revisited  in  a  broader  mathematical 
context  whose  implications  are  of  relevance  to  vision. 
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1.  INTRODUCTION 

A  fundamental  problem  in  early  vision  is  that  of  inferring  the  three-dimensional  geometry 
of  visible  surfaces  in  a  scene  from  the  intensity  information  available  in  the  retinal  images. 
It  seems  that  advanced  biological  vision  systems  subdivide  this  surprisingly  difficult  problem, 
distributing  its  solution  over  a  number  of  subsystems  or  modules  which  contribute  to  the  recovery 
of  information  about  surface  shape.  Each  module  specialises  in  the  interpretation  of  specific  classes 
of  image  cues  and,  to  a  first  approximation,  it  performs  its  task  independently  of  the  other 
subsystems.  Examples  of  primary  modules  which  have  been  identified  include  those  responsible 
for  stereo  vision  and  the  perception  of  motion. 

The  computational  framework  set  forth  by  Marr  [Marr,  1976,  1982;  Marr  and  Poggio,  1977] 
has  had  a  strong  influence  in  the  understanding  of  early  vision.  In  particular,  it  has  provided  a 
paradigm  which  dictates  that  the  modules  initially  be  characterized  in  terms  of  the  visual  tasks 
which  they  must  perform  and,  subsequently,  that  they  be  studied  in  terms  of  the  computational 
processes  through  which  they  perform  these  tasks.  The  various  computational  processes  in  early 
vision  transform  symbolic  representations  of  images  into  symbolic  representations  of  surfaces 
over  several  stages  of  analysis.  An  understanding  of  some  of  the  details  of  these  processes  and 
representations  has  evolved  recently,  especially  at  the  stages  closer  to  the  image,  where  much 
of  the  work  has  been  inspired  by  recent  advances  in  neuroscience  ([Marr,  1976],  [Marr  and 
Hildreth,  1980],  [Marr  and  Poggio,  1979],  [Ullman  1979a],  etc.).  On  the  other  hand,  our  insights 
at  stages  closer  to  explicit  surface  and  volumetric  representations  are  meager.  For  example, 
these  fundamental  questions  remain  open:  first,  how  is  the  information  generated  by  the  various 
modules,  amalgamated  into  representations  of  surfaces  (see,  e.g.,  [Nishihara,  1981])  and,  second, 
how  do  such  representations  give  rise  in  turn  to  representations  of  the  three-dimensional  properties 
of  objects  in  the  scene  (see,  e.g.,  [Brooks,  1981],  [Marr  and  Nishihara,  1978],  and  [Nevatia  and 
Binford,  1977]). 

The  work  described  in  this  paper  is  part  of  ongoing  research  into  the  problem  of  obtaining 
surface  representations  which  will  be  of  use  to  later  processing  stages  in  vision.  In  the  context  of 
Marr’s  framework,  our  goal  is  to  analyze  the  process  through  which  the  sparse  information  retrieved 
by,  say,  stereopsis  or  analysis  of  motion  is  combined  and  transformed  into  full,  retinocentric 
descriptions  of  surface  shape  consistent  with  our  perception  when  we  look  around  us.  In  particular, 
we  will  argue  that  multiple,  full  surface  representations  spanning  a  range  of  resolutions  are  desirable 
and,  indeed,  show  that  they  may  be  generated  as  an  integral  part  of  a  highly-efficient,  multi-level 
surface  reconstruction  algorithm.  Moreover,  our  approach  seems  sufficiently  general  to  allow  several 
classes  of  surface  shape  information  (such  as  local  depth  or  orientation)  provided  by  a  number  of 
vision  modules  to  be  merged  in  a  meaningful  way. 
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1.1.  Motivation  of  the  Multi-Level  Approach 

To  clarify  our  intentions,  we  will  first  examine  Marr’s  framework  for  early  vision  in  some 
detail.  The  framework  is  characterized  by  at  least  three  major  processing  stages,  each  of  which 
transforms  one  retinocentric  representation  into  another  with  the  purpose  of  inferring,  and  making 
explicit,  relevant  information  about  the  surfaces  in  a  scene.  The  first  stage  transforms  the  intensity 
representations  (or  retinal  images)  into  a  primary  representation,  called  the  primal  sketch  [Marr, 
1976].  Changes  in  the  physical  properties  of  surfaces  almost  always  give  rise  to  intensity  changes 
in  the  images,  and  it  is  at  the  level  of  the  primal  sketch  that  the  locations  of  these  changes  are 
made  explicit.  In  the  second  processing  stage,  specialized  processes,  such  as  those  concerned  with 
stereo  and  shape  from  motion,  infer  information  about  the  shape  of  surfaces  from  the  contents  of 
the  primal  sketch.  Since  inferences  can  typically  be  made  only  at  those  locations  which  have  been 
marked  in  the  primal  sketch,  the  information  generated  is  sparse,  and  it  is  collected  into  sparse 
representations  of  surface  shape  that  are  referred  to  as  the  raw  2\-D  sketch.  The  final  stage  is 
one  of  full  surface  reconstruction  in  which  the  sparse  representations  are  transformed  into  a  full 
2 \-D  sketch  containing  explicit  information  about  surface  shape  at  all  points  in  the  scene. 

The  goal  of  the  first  processing  stage  is  the  detection  of  intensity  changes  in  the  image. 
Recently,  Marr  and  Hildreth  [Marr  and  Hildreth,  1980;  Hildreth,  1980]  proposed  a  theory  of  edge 
detection  which  was  inspired  by  existing  neurophysiological  evidence  and  certain  mathematical 
issues.  An  important  aspect  of  this  theory  is  that  intensity  changes  in  the  images  must  be 
isolated  at  different  scales  of  resolution.  Indeed,  there  is  evidence  that  the  human  visual  system 
detects  intensity  changes  over  a  range  of  resolutions  through  the  use  of  up  to  five  independent, 
spatial-frequency-tuned,  bandpass  channels  [Campbell  and  Robson,  1968;  Wilson  and  Giese,  1977; 
Wilson  and  Bergen,  1979;  Marr,  Poggio,  and  Hildreth,  1980].  The  existence  of  these  independent 
primal  sketch  representations  is  a  crucial  factor  which  contributes  to  the  success  of  some  later 
computations  such  as  stereopsis,  as  modeled  by  the  Marr-Poggio  theory  of  stereo  vision  [Marr  and 
Poggio,  1979]  (see  also  [Mayhew  and  Frisby,  1980,  1981]).  According  to  this  model,  the  bandpass 
nature  of  the  channels  leads  to  an  almost  trivial  solution  to  the  stereo  correspondence  problem 
within  the  disparity  range  of  each  channel.  Detailed  depth  information  over  a  wide  disparity  range 
is  obtained  through  a  process  by  which  the  coarser  channels  control  vergence  eye  movements 
that  bring  the  finer  channels  into  alignment  (general  studies  of  vergence  eye  movements  include 
[Riggs  and  Niehl,  1960]  and  [Rashbass  and  Westheimer,  1961]).  On  the  other  hand,  computations 
such  as  motion  correspondence  [Ullman,  1979a],  whose  function  may  not  depend  critically  on 
the  existence  of  multiple  representations,  may  nevertheless  be  operative  at  each  of  the  levels.  It 
seems  likely  in  any  case  that  multiple  sparse  representations  of  surface  shape  that  span  a  range 
of  resolutions  are  generated  by  roost  of  these  modules. 

In  the  context  of  stereopsis,  Grimson  [1981a,  1982a]  pioneered  the  mathematical  theory 
of  the  subsequent  visual  surface  reconstruction  process  which  transforms  the  sparse  surface 
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descriptions  into  full  ones.  He  proposed  that,  before  reconstruction  begins,  the  multiple,  sparse 
depth  representations  output  through  the  different  bandpass  channels  be  combined  into  a  single 
raw  2J-D  sketch  in  a  way  which  maintains  consistency  across  all  scales.  The  raw  2$-D  sketch  then 
contains  sparse  depth  information  at  the  finest  resolution  possible.  Next,  a  single  reconstruction 
process  operating  at  this  finest  level  generates  a  unique  full  2  J-D  sketch  representing  depth 
information  at  high  resolution.  The  steps  are  illustrated  in  Figure  1,  in  which  only  three  bandpass 
channels  are  shown  for  simplicity. 
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A  single  full  surface  representation  at  the  finest  resolution  possible  certainly  captures  all 
of  the  information  provided  by  the  stereo  module  and  it  should,  in  principle,  be  sufficient  input 
to  later  tasks.  Unfortunately,  a  number  of  technical  problems  arise  with  this  simple  approach. 
First,  in  collapsing  the  multiple  depth  representations  into  one  raw  2^-D  sketch,  information 
potentially  useful  in  later  processing  stages  which  are  concerned  with  object- centered  surface 
descriptions,  3-D  models  of  objects,  and  object  recognition  has  been  discarded  prematurely.  It 
now  seems  likely  that  in  order  for  some  of  these  later  stages  to  succeed  and  work  efficiently, 
surface  representations  at  multiple  scales  will  be  necessary,  just  as  they  are  necessary  at  earlier 
stages  such  as  stereopsis.  In  accordance  with  Marr’s  principle  of  least  commitment  [Marr,  1976],  it 
would  be  wasteful  to  discard  information,  prior  to  surface  reconstruction,  which  may  have  to  be 
regenerated  later.  A  second,  and  more  immediately  serious  problem  is  a  consequence  of  the  great 
bulk  of  incoming  information  within  the  large  raw  2^-D  sketch  which  must  be  processed  at  the 
finest  resolution.  Biologically  feasible  surface  reconstruction  algorithms  such  as  those  developed 
by  Grimson  are  extremely  inefficient  at  generating  full  surface  descriptions  when  faced  with  such 
large  representations.  Roughly  speaking,  the  primary  reason  for  this  inefficiency  is  due  to  the 
local  nature  of  the  algorithms  in  question. 

The  above  problems  may  be  avoided  if  the  sparse  representations  are  not  collapsed  into 
a  single  fine  representation.  Instead,  multiple  full  surface  representations  spanning  a  range  of 
resolutions  ought  to  be  generated  by  the  reconstruction  process  itself  and  made  available  to 
processing  stages  beyond.  The  multi-level  surface  reconstruction  algorithm  which  we  will  develop 
in  this  paper  accomplishes  precisely  this.  Because  the  algorithm  exploits  information  available  at 
coarser  resolutions,  its  speed  efficiency  is  dramatically  superior  to  that  of  single  level  reconstruction 
schemes.  Order-of-magnitude  improvements  are  typically  observed  for  surfaces  reconstructed  from 
information  provided  by  stereopsis.  On  the  other  hand,  the  expense  in  space  in  maintaining  all 
the  coarser  representations  is  very  worthwhile  since  it  turns  out  to  be  only  a  fraction  of  that 
required  to  maintain  the  finest  one. 

Figure  2  illustrates  the  multi-level  surface  reconstruction  scheme  and  its  incorporation  into 
stereopsis.  A  fundamental  point  to  realise  about  the  multi-level  approach  in  general  is  that 
information  about  surface  depth,  or  for  that  matter  surface  orientation,  is  provided  in  each  of 
the  channels  (i.e.,  sparse  representations)  by  the  various  vision  modules  and,  as  will  be  shown, 
contributes  in  an  optimal  way  to  the  generation  of  the  hierarchy  of  full  surface  representations. 
The  multi-level  scheme  involves  both  intra-level  processes  which  propagate  information  across  a 
single  representation,  as  well  as  inter-level  processes  which  communicate  between  representations. 
The  inter-level  processes  are  further  classified  into  those  which  transfer  information  from  coarser 
levels  to  finer  ones,  and  those  which  transfer  information  from  finer  levels  to  coarser  ones.  At 
this  point,  we  emphasise  that  multiple  representations  of  consistent  accuracy  can  be  achieved 
only  if  such  a  bi-directional  flow  of  information  is  allowed  to  take  place  between  the  levels.  This 
statement  will  be  substantiated  rigorously  in  a  later  section. 
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Figure  2.  Multi-level  approach  to  »urfaee  reconstruction. 
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If  the  processes  and  representations  envisioned  are  to  be  considered  as  models  of  the 
human  visual  system,  their  form  is  constrained  from  below  by  what  can  be  implemented  in 
neuronal  hardware.  Although  our  incomplete  knowledge  renders  premature  the  formulation  of 
precise  arguments  along  these  lines,  some  constraints  which  presently  seem  compelling,  such 
as  parallelism,  locality  and  simplicity  of  computation,  efficiency,  uniformity,  and  extensibility 
[UUraan,  1979b)  have  been  factors  in  the  theoretical  analysis  of  the  surface  reconstruction  problem 
and  in  the  search  for  algorithms  (similar  constraints  are  issues  when  considering  implementations 
within  parallel,  pipelined  computer  architectures).  Once  a  specific  algorithm  is  selected  based 
on  these  constraints,  it  may  be  implemented  and  its  performance  can  be  evaluated  empirically 
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in  terns  of  the  original  computational  goal.  If  the  performance  is  consistent  with  psychological 
evidence,  the  algorithm  may  be  thought  of  as  constituting  an  existence  proof  that  the  theory 
adequately  models  the  human  vision  system.  This  is  what  we  are  ultimately  striving  for  in  our 
study  of  the  surface  reconstruction  problem. 

1.2.  Overview 

In  this  paper,  we  lay  down  the  mathematical  foundations  of  a  multi-level  approach  to  visual 
surface  reconstruction,  primarily  in  the  context  of  stereo  vision.  With  the  help  of  a  physical 
model,  the  basic  surface  approximation  problem  is  given  an  intuitive  interpretation.  Although  it  is 
in  general  a  nontrivial  matter  to  solve  this  problem,  our  model  suggests  the  application  of  potent 
methods  which  have  arisen  out  of  classical  mathematical  physics  —  the  calculus  of  variations, 
optimal  approximation  theory,  and  functional  analysis.  Aspects  of  the  above  formalisms  are 
employed  to  render  our  problem  amenable  to  solution  by  numerical  techniques.  The  development 
is  as  follows. 

•  In  Chapter  2,  visual  surface  reconstruction  is  cast  as  an  optimal  approximation  problem 
which  involves  finding  the  equilibrium  position  of  a  thin  flexible  plate  undergoing  bending. 
The  problem  is  posed  formally  as  a  variational  principle  which  we  propose  to  solve  by  first 
converting  it  to  discrete  form  using  the  finite  element  method. 

•  In  Chapter  3,  we  prepare  the  way  for  applying  this  discrete  approximation  method  by 
finding  a  set  of  minimal  conditions  for  our  continuous  problem  to  have  a  unique  solution. 
We  show  that  these  requirements  will  almost  always  be  satisfied  in  practice,  so  that  we  can 
consider  our  surface  approximation  problem  to  be  well-posed,  and  can  proceed  to  obtain  the 
solution. 

•  In  Chapter  4,  we  turn  to  the  task  of  converting  our  continuous  problem  into  discrete 
form.  To  do  so,  we  define  a  simple  nonconforming  finite  element  which  will  constitute  the 
basis  of  our  local,  piecewise  continuous  representation  of  surfaces.  Because  the  element  is 
nonconforming,  we  first  must  prove  that  it  leads  to  unique  discrete  approximations,  and 
that  these  approximations  converge  to  the  exact  solution  as  the  elements  decrease  in  size. 
Having  done  this,  we  derive  the  discrete  surface  approximation  problem  as  a  large  system 
of  linear  equations. 

•  In  Chapter  5,  we  face  the  task  of  solving  this  linear  system  efficiently  in  a  biologically-feasible 
way,  and  it  is  here  that  we  motivate  the  multi-level  approach  for  obtaining  the  solution. 
The  approach  involves  setting  up  a  hierarchy  of  discrete  surface  approximation  problems 
which  span  a  range  of  resolutions  and  exploit  the  information  available  at  each  scale,  and 
subsequently  invoking  a  multi-level  algorithm  to  solve  them  simultaneously.  We  demonstrate 
the  efficient  performance  of  the  multi-level  surface  reconstruction  algorithm  on  constraints 
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from  stereopsis,  and  demonstrate  that  it  generates  a  useful  hierarchy  of  accurate  surface 
representations. 

•  In  Chapter  6,  we  reexamine  our  surface  reconstruction  problem  and  show  that  it  is  a  special 
case  within  a  general  class  of  optimal  interpolation  problems  involving  arbitrary  degrees  of 
continuity,  in  any  number  of  dimensions.  These  general  problems  involve  the  minimisation 
of  functionals  which  possess  a  number  of  invariance  properties  making  them  attractive  for 
application  to  problems  in  early  vision  whose  solutions  require  the  iterative  propagation  of 
smoothness  constraints  across  retinocentric  representations. 

•  In  Chapter  7,  we  conclude  by  discussing  the  overall  implications  of  our  approach  to  issues 
concerning  the  isolation  of  depth  discontinuities,  and  the  incorporation  of  other  sources  of 
information  such  as  surface  orientation.  We  discuss  possible  solutions  to  these  problems  in 
view  of  our  finite  element  representation  of  surfaces  and  the  multi-level  surface  reconstruction 
algorithm. 

•  For  convenience,  in  two  of  the  appendices,  we  cover  the  relevant  mathematical  background 
of  the  finite  element  method  and  the  iterative  solution  of  large  linear  systems. 

2.  THE  MOST  CONSISTENT  SURFACE 

The  sparse  information  about  surface  shape  retrieved  by  the  various  vision  modules  is  in 
general  underconstraining.  That  is  to  say,  it  is  insufficient  to  compel  a  unique  inference  of  the 
physical  properties  of  the  surfaces  in  the  scene.  Yet,  even  when  presented  with  impoverished 
stimuli  (such  as  random  dot  stereograms  [Julesz,  1971]  or  kinetic  depth  effect  displays  [Wallach 
and  O’Connell,  1953;  Wallach,  1959;  Johansson,  1975;  Ullman,  1979a;  etc.]),  the  human  visual 
system  routinely  arrives  at  unique  interpretations  and  our  typical  perception  is  a  stable  one  of 
full  surfaces  in  depth.  Clearly,  the  visual  system  must  invoke  certain  assumptions  which  provide 
enough  additional  constraint  to  allow  a  unique  full  surface  representation  to  be  computed  from 
the  sparse  information  provided.  However,  these  additional  assumptions  must  be  plausible  in  that 
they  reflect  certain  domain-dependent  expectations.  For  example,  in  stereo  vision,  the  sparse 
information  takes  the  form  of  depth  constraints  which  embody  measurements  of  the  distances 
from  the  viewer  to  the  surfaces  of  objects  in  the  scene.  The  additional  assumptions  should  then 
be  based  on  general  expectations  about  physical  properties  of  surfaces  in  the  visual  world,  as  well 
as  aspects  of  the  optical  and  computational  processes  taking  part  in  the  generation  of  the  depth 
constraints. 

Grimson  [1981a]  explored  a  number  of  issues  along  these  lines.  Qualitatively,  his  thesis  is 
as  follows.  A  surface  in  the  scene  which  varies  radically  in  shape  usually  gives  rise  to  intensity 
changes  which  are  marked  in  the  primal  sketches  as  zero-crossings  of  the  Laplacian  of  the 
Gaussian-convolved  images  (VSG  *  I  —  see  [Marr  and  Hildreth,  1980;  Hildreth,  1980]).  Moreover, 
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it  it  only  at  the  locations  of  aero-crossings  that  the  Marr-Poggio  stereo  algorithm  can  generate 
measurements  of  the  distance  to  the  surface,  in  the  form  of  explicit  depth  constraints.  Therefore, 
the  surface  cannot  in  general  be  varying  radically  in  depth  between  the  constraints  to  which  it 
gave  rise.  By  introducing  this  additional  surface  smoothness  assumption,  the  goal  of  accurately 
reconstructing  the  shape  of  visible  surfaces  and  thereby  computing  full  surface  representations 
consistent  with  our  perception  is  attainable  in  principle.  A  theoretical  proof  of  this  statement  lies 
in  the  domain  of  mathematics.  In  the  next  section,  we  take  the  first  step  by  rigorously  formulating 
the  surface  reconstruction  problem  as  :in  optimal  approximation  problem  in  which  the  smoothness 
assumption  has  a  clear  intuitive  interpretation,  and  which  eventually  leads  to  efficient  multi-level 
algorithms  for  its  solution. 

2.1.  A  Physical  Interpretation  —  The  Bending  of  a  Thin  Plate 

Visual  surface  reconstruction  can  be  characterized  formally  aB  a  constrained,  optimal 
approximation  problem  in  two  dimensions.  In  the  context  of  stereo  vision,  where  constraints 
embody  depth  measurements  to  surfaces  in  the  scene,  the  goal  is  to  reconstruct,  as  accurately  as 
possible,  the  shape  of  the  surface  which  gave  rise  to  these  measurements.  Of  course,  it  is  necessary 
that  we  be  able  to  deal  with  the  complication  arising  from  arbitrarily-placed  constraints  since,  as 
has  been  argued,  constraints  of  this  type  are  generated  naturally  by  stereopsiB  and  other  vision 
modules.  More  rigorously,  the  problem  can  be  stated  as  follows:  given  a  finite  set  of  arbitrarily 
located  distinct  point  constraints  on  a  plane,  each  constraint  having  a  real  scalar  value  associated 
with  it,  find  the  unique  optimal  function  of  two  variables  which  is  most  consistent  with  these 
constraints.  Our  notion  of  consistency  will  be  defined  shortly.  We  consider  the  solution  to  our 
problem  to  be  a  full  surface  representation  in  that  it  makes  explicit  our  best  estimate  of  the 
distance  to  every  visible  point  on  the  surface  in  the  scene. 

The  constraints  provided  by  the  stereo  computation  are  never  completely  reliable.  Errors  due 
to  noise,  and  errors  in  matching  corresponding  zero-crossings  are  bound  to  occur.  This  suggests 
that  we  should  not  try  to  interpolate  the  given  data  exactly  because  a  few  “bad”  constraints 
can  have  a  detrimental  effect  on  the  shape  of  the  recovered  surface.  Relaxing  the  interpolation 
requirement  turns  our  problem  into  one  of  surface  approximation  in  which  we  would  like  to 
maintain  control  over  how  closely  the  surface  fits  the  data. 

By  thinking  in  terms  of  an  optimal  surface,  we  imply  the  choice  of  a  suitable  criterion  that 
will  allow  us  to  measure  the  optimality  of  admissible  functions.  A  suitable  criterion  for  measuring 
the  optimality  of  surfaces  in  the  context  of  surface  approximation  in  stereopsis  translates  into 
a  precise  mathematical  statement  which  captures  intuitive  notions  about  the  smoothness  of 
admissible  surfaces  as  well  as  their  closeness  of  fit  to  the  known  depth  constraints.  Perhaps  the 
intuitively  clearest  treatment  of  our  problem  is  in  terms  of  a  physical  model.  Consider  a  planar 
region  U,  the  region  within  which  we  wish  to  obtain  an  optimal  approximating  surface  most 
consistent  with  a  finite  set  of  sparse  constraints.  Let  ub  imagine  that  the  constraints  constitute  a 
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Figure  3.  The  physical  model  for  surface  approximation. 


set  of  vertical  pins  scattered  inside  Cl,  the  height  of  an  individual  pin  being  related  to  the  distance 
from  the  viewer  to  the  surface  in  the  scene.  Suppose  that  we  take  a  thin  flexible  plate  of  clastic 
material  that  is  planar  in  the  absence  of  external  forces,  and  constrain  it  to  pass  near  the  tips  of 
the  pins  by  attaching  ideal  springs  between  the  pin  tips  and  the  surface  of  the  plate  as  shown  in 
Figure  3.  It  is  not  difficult  to  imagine  the  equilibrium  position  of  the  plate  as  a  function  of  the 
various  pin  heights. 

Intuitively,  the  equilibrium  position  of  the  thin  plate  is  a  “fair”  approximating  surface  in 
that  it  will  exhibit  a  sufficient  amount  of  smoothness  between  the  constraints.  Moreover,  on 
quantitative  grounds,  there  is  evidence  to  suggest  that  such  a  surface  is  indeed  an  optimal  one  in 
terras  of  the  imaging  process  [Grimson,  1982b].  In  any  case,  we  have  a  reasonable  physical  model 
for  the  optimal  approximating  surfaces  and,  moreover,  this  model  will  suggest  good  strategies  for 
solving  our  problem. 

We  emphasize  however  that  the  appropriateness  or  the  model  depends  on  two  important 
issues.  The  first  involves  ensuring  that  a  unique  solution  exists,  and  the  second  is  to  guarantee 
that  the  solution  is  meaningful  in  view  of  the  constraints.  Firstly,  we  realize  that  the  plate-spring 
system  will  be  unstable  for  certain  pin  configurations.  If  we  have  but  a  single  pin,  then  a  stable 
equilibrium  does  not  exist,  as  the  plate  has  two  unconstrained  degrees  of  freedom  (rotation  about 
the  axis  of  the  pin  is  excluded).  A  similar  degenerate  situation  arises  for  the  case  of  any  number 
of  pins  arranged  linearly,  the  plate  then  having  one  unconstrained  degree  of  freedom.  Clearly 
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at  least  three  aoncollinearly  arranged  pins  are  required  to  assure  that  a  unique  state  of  stable 
equilibrium  exists.  Secondly,  a  reasonable  choice  must  be  made  for  the  stiffness  of  the  springs. 
If  the  springs  are  too  strong  in  relation  to  the  rigidity  of  the  plate  material,  then  a  pin  whose 
height  deviates  significantly  from  that  of  its  neighbors  (i.e.  an  erroneous  constraint)  will  place  an 
abnormally  large  amount  of  strain  on  the  plate  locally  and  have  undesirable  effects  on  the  shape 
of  the  surface.  On  the  other  hand,  if  the  springs  are  too  weak,  the  intrinsic  rigidity  of  the  plate 
can  overwhelm  them  and  the  plate  will  remain  nearly  planar  over  large  variations  in  the  height 
of  the  pins.  In  the  limit  of  a  rigid  plate,  the  resulting  planar,  least-squares  approximation  would 
be  meaningless  in  that  the  solution  does  not  in  general  lie  close  to  constraints  other  than  those 
arising  from  nearly  flat  surfaces. 

We  will  now  proceed  to  a  mathematical  characterization  of  the  above  physical  model.  To  do 
so,  we  apply  the  well-known  minimum  potential  energy  principle  from  classical  mechanics,  which 
states  that  the  potential  energy  of  a  physical  system  in  a  state  of  stable  equilibrium  is  at  a  local 
minimum.  For  the  model,  the  potential  energy  in  question  is  that  due  to  deformation  of  the  plate 
and  springs,  as  well  as  the  energy  imparted  by  any  externally  applied  forces.  In  this  sense  then, 
the  surface  we  seek  is  one  of  minimal  energy. 

First,  consider  the  plate.  It  is  known  (see,  e.g.,  [Courant  and  Hilbert  1953,  pg.  250],  [Landau 
and  Lifihitz,  1970])  that  the  potential  energy  of  a  thin  plate  under  deformation  is  given  by  an 
integral  of  a  quadratic  form  in  the  principle  curvatures  of  the  plate.  If  the  principle  curvatures 
of  the  deformed  plate  are  denoted  by  tc\  and  tc-i,  the  potential  energy  density  is  given  by  an 
expression  of  the  form 

7j- (*?  *2)  ~f  Bk.\ K2  =  2A ^  1  — ( A — B)k\K2, 

where  A  and  B  are  constants  determined  by  the  plate  material.  The  expression  J(«i  4-  *2)  is  the 
first  or  mean  curvature  and  is  the  second  or  Gaussian  curvature  of  the  plate’s  surface  (see, 
e.g.,  [Hilbert  and  Cohn-Vossen,  1952]). 

Let  the  function  v(x,  y)  denote  the  deflection  of  the  plate  normal  to  the  region  fl  which  can 
be  taken  to  lie  in  the  X-Y  plane.  Assuming  that  the  deflection  function  and  its  partial  derivatives, 
v,vx,vv, . . .  are  small,  it  can  be  shown  (Bee  e.g.  [Rektorys,  1969,  pg.  368])  that 

*1  +  *2  1  A  2 

- -  »  -  Aw,  K.  j  K 2  sa  VxzVvy  —  V‘v, 

where  A  =  ^  denotes  the  Laplacian  operator.  Thus,  the  potential  energy  density  can  be 

written  in  the  following  forms: 
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e,  =  -(At/)2  -  (1  -  <r)[vxxvvy  -  ua„) 


\(VIZ  H  at/,,,  +  V^)  +  a(VItVyy  -  V\y) 

(Av)a-f(l-a)(t,L  +  2<+vJy)]I 


(1) 


apart  from  a  multiplicative  constant  which  depends  on  the  physical  properties  of  the  elastic 
material  of  the  plate,  and  which  has  been  set  to  unity  without  loss  of  generality.  The  constant  a, 
called  the  Poisson  ratio,  measures  the  change  in  width  as  the  material  is  stretched  lengthwise.1 
The  desired  potential  energy  of  deformation  is  obtained  by  integrating  the  energy  density  over 
the  domain  in  question,  and  is  given  by 

=  J  fn\(Av)3  ~  0  ~  ff)(vzzv»v  ~  viy)  dxdy. 

To  the  potential  energy  of  deformation,  we  must  add  the  potential  energy  due  to  any 
external  forces  which  may  be  present.  The  energy  due  to  a  force  density  g( x,y)  applied  to  the 
surface  of  the  plate  (such  as  the  effect  of  gravity)  is  given  by 

£2(u)  =  —  j  J  gvdx  dy. 

External  forces  and  bending  moments  may  also  be  applied  around  the  boundary  3H  of  0.  The 
energy  due  to  a  force  density  p(s)  on  the  boundary  (s  denotes  arc  length  along  the  boundary)  is 

^3  (*>)  =  —  /  P(s)v<k, 

Jan 

while  the  energy  due  to  bending  moments  applied  around  the  boundary  is 

f«(v)  =  -  /  rf*> 

Jan  on 


where  m(s)  is  the  density  of  applied  bending  moments  normal  to  the  curve  and  denotes  the 
directional  derivative  along  the  outward  normal  to  dfl. 


Finally,  we  must  account  for  the  potential  energy  of  deformation  of  the  springs.  Let  C 
denote  the  set  of  points  in  fl  at  which  the  imaginary  pins  are  located;  that  is,  the  sparse  set  of 
locations  at  which  the  surface  is  constrained.  Furthermore,  denote  the  height  of  the  pin  (the  value 
of  the  constraint)  by  real  scalars  C(z,,y.)  and  the  stiffness  of  the  spring  attached  to  it  (influence 
of  the  constraint)  by  positive  constants  P{z,,Vi)  for  all  (x,,  y.)  €  C  .  According  to  Hooke’s  law  for 
an  ideal  spring,  the  total  potential  energy  of  deformation  in  the  springs  is  given  by 


'From  the  last  expression  in  (1),  it  is  apparent  that  the  potential  energy  density  may  be  considered 
to  be  a  convex  combination  with  parameter  o  of  the  square  of  the  Laplacian  ard  a  quadratic  term  in 
second-order  partial  derivatives  of  the  deflection  function,  (vZI  +  1v\t  vy>).  This  fact  will  be  used  in 
a  subsequent  discussion. 
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£s(v) — r  Vt)  c(*i.v>)]  • 

(*i.Vi)€C 

The  equilibrium  state  of  the  mechanical  system  can  be  obtained  as  the  solution  to  the 
following  minimisation  problem  which  is  referred  to  as  a  vartattonal  principle: 

The  deflection  of  the  plate  at  equilibrium  is  that  function  u  from  a  set  V  of  admissible 

functions  v,  for  which  the  total  potgiiiiaLuuxqy - — 

£(v)  —  £,(«)  +  £2(0)  +  £3(v)  +  £e(v)  4-  £>(«)  (2) 

is  minimal. 

Thus,  quantitatively,  the  "most  consistent”  surface  which  we  seek  is  the  one  having  minimal 
energy  f . 

The  visual  surface  approximation  problem  has  been  posed,  in  integral  form,  as  a  variational 
principle  which  is  quadratic  in  that  it  involves  terms  that  are  at  most  quadratic  in  v  and  its 
derivatives.  Through  the  formalism  of  the  calculus  of  variations  one  can  express  the  necessary 
condition  for  the  minimum  as  the  Euler-Lagrange  equation  which  in  thiB  case  can  be  shown  to 
be  a  linear,  self-adjoint,  partial  differential  equation.  Much  of  classical  mathematical  physics  is 
based  on  this  duality,  and  it  provides  numerous  techniques  for  solving  our  problem,  those  which 
are  directed  towards  the  variational  principle,  as  well  as  those  which  are  directed  towards  its 
Euler-Lagrange  equation. 

Whatever  the  strategy,  although  it  is  conceivable  that  the  exact  analytical  solution  u 
could  be  derived,  it  is  normally  impossible  to  do  so  for  all  but  the  Bimplest-shaped  domains  fl. 
Consequently,  we  are  led  to  consider  a  numerical  approach  in  which  we  somehow  convert  our 
continuous  problem  into  a  discrete  problem  whose  numerical  solution  closely  approximates  the 
exact  continuous  solution  u.  We  propose  to  employ  what  is  probably  the  most  potent  tool  for 
obtaining  discrete  approximations  currently  available  —  the  finite  element  method.  The  method  is 
applied  to  the  variational  principle  directly  and,  because  the  variational  principle  is  quadratic,  the 
resulting  discrete  problem  will  take  the  particularly  simple  form  of  a  linear  system  of  algebraic 
equations.  The  main  advantage  of  the  finite  element  method  is  its  generality.  In  the  context 
of  our  surface  approximation  problem,  it  can  be  applied  over  domains  fl  of  complicated  shape, 
and  it  is  not  limited  to  uniform  discretizations  of  these  domains.  The  importance  of  the  latter 
property  in  the  context  of  vision  is  evident  when  one  considers,  for  example,  the  nonuniform 
structure  of  the  retina  where  it  is  known  that  resolution  decreases  approximately  linearly  with 
eccentricity  (see  [Wilson  and  Giese,  1977]  and  [Wilson  and  Bergen,  1979]  for  a  quantitative  model 
of  this  phenomenon  in  terms  of  the  spatial-frequency  channels  in  early  vision).  Moreover,  the  finite 
element  method  leads  to  linear  systems  which  are  readily  solvable  in  a  parallel,  iterative  fashion 
by  a  sparsely-interconnected  network  of  simple  processors,  a  mechanism  which  seems  prevalent  in 
early  vision. 
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For  leveral  reasons,  we  have  avoided  the  alternate  route  of  using  the  well-known  finite 
difference  method  to  discretize  the  associated  Euler-Lagrange  equations  (see,  e.g.,  [Collatz,  1966], 
[Forsythe  and  Wasow,  1960],  [Smith,  1977]).  The  finite  difference  method  is  much  more  restrictive 
in  that  it  practically  limits  us  to  uniform  discretisations,  the  underlying  convergence  theory  it 
much  less  well  developed  and,  perhaps  most  importantly,  it  becomes  very  difficult  to  discretize 
the  natural  boundary  conditions  associated  with  our  surface  approximation  problem,  a  task  which 
is  done  trivially  by  the  finite  element  method. 

The  mathematical  background  of  the  finite  element  method  that  is  of  relevance  to  our 
problem  is  included  in  Appendix  A  for  convenience.  The  appendix  introduces  the  required  theory 
and  lists  the  fundamental  theorems  which  we  will  invoke  in  applying  the  method  to  the  task 
at  hand.  The  process  will  consist  of  several  steps.  First  we  pose  the  variational  principle  in  an 
abstract  form  that  is  the  basis  of  the  mathematical  machinery  presented  in  Appendix  A.  Next, 
we  determine  formally  the  requirements  on  the  boundary  conditions  that  must  be  satisfied  to 
ensure  that  our  variational  principle  is  well-posed;  i.e.,  has  a  unique  solution.  Only  then  can  we 
proceed  to  apply  the  finite  element  method  to  approximate  the  solution. 

3.  THE  VARIATIONAL  PRINCIPLE 


In  this  chapter  we  analyze  the  continuous  variational  principle  which  embodies  our  visual 
surface  approximation  problem,  in  preparation  for  the  application  of  the  finite  element  method. 
In  view  of  the  formalism  presented  in  Appendix  A,  our  first  goal  is  to  state  the  variational 
principle  in  the  abstract  form;  that  is,  to  isolate  the  energy  inner  product  which  characterises  our 
minimization  problem.  We  then  derive  the  associated  Euler-Lagrange  equation  and,  in  the  process, 
consider  the  various  forms  of  boundary  conditions  that  can  be  imposed.  Finally,  we  choose  the 
appropriate  form  of  these  conditions  in  view  of  our  visual  surface  approximation  problem  and 
obtain  formally  the  minimum  requirements  for  our  variational  principle  to  be  well-posed. 

3.1.  The  Energy  Inner  Product 

According  to  equation  (2.2),  our  variational  principle  asserting  that  the  equilibrium  state  of 
the  thin  plate  is  a  minimal-energy  configuration,  may  be  stated  mathematically  as  the  minimisation 
of  the  expression 

£(v)=  f  f  i(Av)2  —  (1  —<r)(vxxvvy  -«*„)  —  gvdxdy 

02  ,  f  dv  8  i  (i) 

EM*"  • 

Here,  we  have  assumed  that  the  spring  stiffnesses  P(xuVi)  =  for  all  (x,,  y<)  €  C .  The  admissible 
space  V  for  our  variational  principle  is  in  general  a  subspace  of  the  second-order  Sobolev  space 
l/2(fi)  over  the  region  H  (refer  to  the  discussion  in  Appendix  A).  If  u  €  V  minimises  C,  then 
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£(u)  <  if(u-f-ev)  for  any  v  £  V.  Therefore,  to  obtain  the  necessary  condition  for  the  minimum,  we 
substitute  for  v  the  small  variation  u  -f-  tv  about  u  and  equate  to  sero  the  derivative  with  respect 
to  e  for  e  =  0.  Equivalently,  we  may  perform  the  variation  using  the  rules  of  differentiation: 


$<f(ti)  =  J  J^AuS(Au)  —  (l—a)[uTZS(uyy)  +  uvv6(Mzz)—2MZyS{u,Zy)\  —  g6udxdy 

—  I  p(a)Suda—  f  m(a)S\  ~  ) da  +  p  V  —  Cf*. *•)]*“(*«* »<)• 

Jan  Jan  Vdn/  , 


Since  variation  commutes  with  differentiation ,  S(Au)  =  A (6u),  &(§%)  =  ^(5u),  S(u yy)  =  (S u)yy, 
etc.  If  we  now  let  v  =  6 u,  and  set  ££  =  0,  we  obtain 


n 


AuAv  —  (1  —  tr)(uzzvvv  -f  riyyVzx  —  2uivuiy)  —  gvdxdy 


~  fan^Vda~fanm^$ndS  +  l3  ~  c(*.v<>M*m Vt)  =  0. 


(2) 


(*s»y<)GC 


Equations  (1)  and  (2)  may  be  cast  in  our  abstract  variational  formulation  of  Appendix  A. 
The  key  is  in  identifying  the  energy  inner  product  as  the  bilinear  form 


a(u,  v)  ==  J  I  AuAv  —  (1  —  <r){uzzVyy  +  UyVvzz  —  2ux„uly)didy  +  p  ^Z  u(x*>  V>)> 

n  (x.,y.)ec 

(3) 

and  in  writing  the  linear  form  as 


•  (4) 


Clearly  then,  (1)  asserts  that  we  are  to  minimise  the  quadratic  functional  S(v)  =  $a(t>,  v)  —  f(v), 
as  required  in  the  definition  of  the  abstract  variational  principle  (Definition  A.l).  On  the  other 
hand,  (2)  which  expresses  the  necessary  condition  for  the  vanishing  of  the  first  variation  may  be 
written  as  a(u,  v)  =  /(v),  as  expected  from  the  discussion  leading  up  to  the  variational  equation 
(A.  10). 

Having  obtained  expressions  for  the  bilinear  and  linear  forms,  we  can  proceed  to  bring  the 
finite  element  method  to  bear  on  the  problem.  Before  doing  so,  however,  it  is  imperative  that 
we  carry  the  analysis  further  so  that  we  can  express  the  necessary  condition  for  a  minimum 
as  a  partial  differential  equation,  explore  the  issue  of  boundary  conditions,  and  ensure  that  the 
problem  is  well-posed. 


3.2.  The  Euler-Lagrange  Equation  and  Boundary  Conditions 

For  the  duration  of  this  section,  we  will  ignore  the  summation  term  arising  from  the  (spring) 
constraints,  since  its  presence  will  complicate  the  notation  while  making  no  significant  contribution 
to  the  discussion.  First,  we  will  transform  the  energy  inner  product  o(-,  -)  given  in  (3)  using 
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integration  by  parts  in  two  dimensions;  i.e.  Green’s  theorem.  Let  n  be  the  outward  normal  vector 
to  3fl,  t  be  the  usual  tangent  vector  along  dfl,  and  ^  and  $  denote  partial  differentiation 
along  the  normal  and  tangent  respectively.  Assuming  that  u  is  fourth-order  differentiable,  Green’s 
identity  (see,  e.g.,  [Ciarlet,  1978,  pg.  14],  [Rektorys,  1969]) 


vAudxdy 


-  f 

Jan  dn 


da 


may  be  used  to  transform  the  term  AuAu  arising  from  the  mean  curvature  of  the  surface: 


J  L  ' dy  =  /  L vAiudx  dy  +  L  Au^  da  -  L  Wli(A,l) *• 


(5) 


where  A2u  =  AAu  =  uUi,  -(-  2uIlyy  -(-  uyyyy.  On  the  other  hand,  the  Gaussian  curvature  term 
can  be  transformed  using  the  identity  [Ciarlet,  1978,  pg.  15] 


//„*■ 


+  ti; 


yyVxx 


2UZyVZy 


***>  =  !  Sr*-/ 

Jon  dt2  dn  Jo 


d2u  dv 
on  dndt  dt 


da. 


(6) 


If  the  boundary  is  sufficiently  smooth,  it  can  be  shown  (see  [RektoryB,  1980,  pp.  268-269])  that 
the  second  boundary  integral  can  be  written  as 


L 


d2u  dv 
on  dndt  dt 


Substituting  equations  (5)-(7)  into  (3)  (and  ignoring  the  constraint  term),  we  obtain 

a(u,  v)  =  J  J^vA2udxdy  +  J  ^P(u)v  da  +  J^M(u)^ds, 


(7) 


where 


M(u)  =  Au-(1— o)|^. 

Thus,  the  necessary  condition  for  the  minimum  (2)  becomes 

J  (A2u  -  g)vdx  dy  -f  [P(u)  —  p(s)]v  da  +  |A/(«)  —  m(a)}~da  =  0 

Now,  since  the  above  equation  must  bold  and  since  v  and  §%  are  arbitrary  on  the  dosed 
region  fi,  we  must  have 


A2u  = 


in  f). 


(8) 


This  is  the  fourth-order,  linear,  self-adjoint  Euler-Lagrange  equation  that  governs  the  small 
deflection  of  a  thin  plate  at  equilibrium,  and  it  is  satisfied  by  u  inside  fl  regardless  of  the 
boundary  conditions  on  dfl.  In  its  homogeneous  form  A2u  =  0,  it  is  called  the  biharmonic 
equation.  Furthermore,  u  must  satisfy  the  natural  boundary  condition* 
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P(u)  —  p(a)  and  M(u)  =  m(a)  on  dft.  (9) 

According  to  (6),  the  integral  over  (1  of  the  Gaussian  curvature  approximation  (vxxvvv  —  2 
has  no  effect  whatsoever  on  the  Buler-Lagrange  equation,  but  contributes  only  to  the  boundary 
conditions.2  This  reflects  the  fact  that  the  Gaussian  curvature  of  a  surface  patch  is  invariant  as  the 
surface  is  subjected  to  arbitrary  bending,  and  that  its  average  value  over  the  patch  depends  only 
on  the  tangent  planes  of  the  surface  around  the  periphery  of  the  patch  [Hilbert  and  Cohn-Vossen, 
1952,  pp.  193-204],  This  invariance  property  of  the  Gaussian  curvature  renders  it  inappropriate  as 
a  measure  of  surface  consistency.  For  example,  it  cannot  distinguish  two  developable  surfaces  such 
as  a  wildly  undulating  sinusoidal  surface  and  a  planar  surface,  both  of  which  are  cylinders  and 
therefore  have  sero  Gaussian  curvature  everywhere.  On  the  other  hand,  the  mean  curvature  does 
not  in  general  remain  invariant  under  bending  and  therefore  plays  a  vital  role  in  our  energy  inner 
product.  This  is  evident  from  equation  (1)  —  no  value  of  a  can  make  (At/)2,  which  approximates 
the  mean  curvature,  vanish.8 

Another  consequence  of  the  necessary  condition  for  the  minimum  is  that  the  form  of  the 
natural  boundary  conditions  satisfied  by  u  are  determined  by  any  essential  boundary  conditions 
which  may  be  imposed  on  v.  In  general,  we  can  impose  up  to  two  essential  boundary  conditions, 
one  on  o  and  the  other  on  First,  consider  the  case  of  a  simply  supported  plate  where  the 
essential  boundary  condition  ti  =  0  is  imposed  on  dil  but  is  left  unconstrained.  The  solution 
u  must  then  still  satisfy  the  second  condition  in  (9).  We  therefore  have  the  Neumann  boundary 
conditions 

u  =  0,  Af(u)  =  m(s)  on  flfl, 
and,  moreover,  the  first  contour  integral  in  (1)  vanishes. 

Next,  suppose  that  we  also  set  =  0  on  dtl.  Then, 

du 

u  =  —  =  o  on  an,  (io) 

which  are  the  Dirichlet  boundary  conditions  for  the  clamped  plate.  In  this  case,  both  contour 
integrals  in  (1)  vanish  and,  moreover,  a  is  arbitrary  since  it  does  not  appear  in  the  Euler- Lagrange 
equation  (8),  but  only  in  the  natural  boundary  conditions  (9)  which  have  now  been  replaced 
by  (10).  We  can  therefore  greatly  simplify  the  variational  integral.  In  particular,  the  functional 
minimization  problems  involving 


‘Expressions  possessing  this  property  are  called  divergence  expressions  [Courant  and  Hilbert,  1953]. 

‘Brady  and  Horn  [1981,  pg.  29]  state  that  “the  choice  of  which  performance  index  to  use  is  reduced 
to  the  square  Laplacian,  the  quadratic  variation,  and  linear  combinations  of  them”.  We  stress  that  one 
should  be  careful  not  to  choose  that  linear  combination  which  results  in  a  divergence  expression  (the 
Gaussian  curvature)  and  therefore  hr-i  an  identically  sero  Euler-Lagrange  equation.  Recall  from  equation 
(2.1)  that  the  small-deflection  theory  of  the  thin  plate  allows  only  a  convex  combination  so,  fortunately, 
it  is  free  from  danger,  in  this  respect. 
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S(v)  =  J  ^(Av)a  —  gvdxdy,  for  a=  1, 

and 

£(«)  =  /  /;  +  2w*u  +  vly)  ~  gvdxdy,  for  o  =  0, 

are  equivalent  in  the  Dirichlet  case. 

Finally,  consider  the  case  of  a  free  boundary;  that  is,  when  the  externally  imposed  force  p(s) 
and  moment  m(s)  on  the  boundary  are  zero.  Then,  there  are  no  constraints  on  v  but,  according 
to  (9),  u  must  satisfy 

P(u)  =  M (u)  =  0  on  SO. 

These  are  the  natural  boundary  conditions  satisfied  by  the  solution  for  the  case  of  the  free  plate. 
Once  again,  the  contour  integrals  in  (1)  vanish  and  the  energy  functional  takes  the  simple  form 

f  H  =  /  ja  ^(At,)S  —  0  —  -  v2xv )  —  gv  dx  dy. 

In  general,  the  admissible  space  V  is  the  subspace  of  the  Sobolev  space  U2((l)  which  satisfies 
the  essential  boundary  conditions  imposed  on  the  plate.  If,  for  example,  a  portion  of  the  edge 
of  the  plate  is  simply  supported,  V  will  consist  of  functions  which  satisfy  the  essential  condition 
v  —  0  on  that  portion  of  dfl.  If  part  of  the  edge  of  the  plate  is  damped,  then  v  =  =  0  on 

that  part  of  dfl.  On  the  other  hand,  if  part  of  the  edge  is  free,  then  no  constraints  are  imposed 
on  v  over  that  portion  of  the  boundary  and,  in  the  case  of  the  free  plate,  V  =  J/2(fl).  Of  course, 
the  plate  cannot  be  “too  free”  on  fi,  because  then  the  physical  system  cannot  achieve  stable 
equilibrium  and  a  unique  solution  would  not  exist.  Precisely  how  much  freedom  can  be  allowed 
will  be  established  formally  in  the  next  section. 

3.3.  When  is  the  Problem  Well-Posed? 

Turning  to  our  visual  surface  approximation  problem,  we  should  at  this  point  choose  the 
appropriate  form  of  boundary  conditions  on  dfl.  Since  the  only  information  about  the  surface 
that  is  provided  by  the  stereo  module,  for  example,  is  embodied  in  the  sparse  constraints,  the 
strategy  of  least  commitment  is  to  “assume  nothing”  about  the  boundaries  of  the  surface.  In 
terms  of  our  plate  problem,  this  means  that  we  should  impose  no  essential  boundary  conditions 
on  the  plate;  that  is,  we  solve  the  fret  plate  problem  whose  admissible  space  V  =  li2(fi). 

If  the  boundary  of  the  plate  is  free,  it  is  clear  that  the  constraints  will  play  a  critical  role  in 
providing  a  unique  state  of  stable  equilibrium.  Our  goal  in  this  section  is  to  specify  the  existence 
and  uniqueness  requirements  mathematically  as  conditions  under  which  the  surface  approximation 
problem  is  well-posed.  To  do  this,  we  will  invoke  Theorem  A.l,  and  satisfy  its  conditions  by 
proving  the  following  two  propositions. 
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Propoiition  1.  The  energy  inner  product  o(-,  •)  is  symmetric. 

Proof.  o(u,  v)  =  a(v,  u)  is  evident  by  inspection  of  equation  (3).  | 

Proposition  2.  If  the  set  of  constraints  C  contains  (at  least)  three  noncolline&r  points,  then  o(-,  ■) 
is  V-elliptic  for  0  <  cr  <  1. 

Proof.  We  want  to  show  that  there  exists  an  a  >  0  such  that  a(v,  v)  >  a||v||s,  for  all  w  £  V. 
To  do  so,  it  is  sufficient  to  show  that  a(v,  v)  =  0  only  if  v  =  0.  We  rewrite  a(v,  v)  as 

o(w,v)  =  /  f  o(Au)2  +(1  —  »)(«?,  +  2v*y  -f  v\y)dxdy  +  0  «(*.,  y.)2- 

°  (*..V.)6  C 

Now,  Av  =  0  only  if  v  is  a  harmonic  function,  while  -f-  2v*y  -f  t>yy)  =  0  only  if  v  is  a  first 
degree  polynomial  (as  can  easily  be  shown  by  integration),  which  is  a  subclass  of  the  harmonic 
functions.  Thus,  the  integral  is  >  0  for  0  <  <T  <  1  and  it  is  sero  only  if  v  is  a  linear  function 
over  fl.  On  the  other  hand,  since  /3  is  poai**?e  by  definition  the  sum  is  also  >  0  and  it  is  sero 
only  if  v(xi,y,)  =  0  for  all  (x,,  y,)  6  C.  Therefore,  if  C  contains  three  noncollinear  points,  then 
o(v,  ti)  =  0  only  if  v  =  0,  implying  that  ^  ,  •)  is  -eiiiptic.  | 

By  Propositions  1  and  2  and  Theorei-  v.  we  are  assured  that  the  continuous  approximation 
problem  is  well-posed  if  0  <  a  <  1  «.ad  the  set  of  constraints  includes  at  least  three  noncollinear 
points.  The  condition  on  the  constraintr  is  not  unexpected  in  view  of  the  arguments  made  in 
Section  2.1.  Physically  speaking,  all  unconstrained  degrees  of  freedom  of  the  plate  must  be 
precluded,  and  three  noncollinear  constraints  is  clearly  the  minimum  requirement  for  this  to  be  the 
case.  In  application  to  natural  images,  the  stereo  algorithm  will  almost  always  generate  at  least 
three  noncollinear  points,  so  we  can,  for  all  practical  purposes,  consider  our  surface  approximation 
problem  to  be  well-posed  so  long  as  0  <  cr  <  1. 

4.  OBTAINING  THE  DISCRETE  PROBLEM 

So  far,  we  have  been  dealing  with  the  continuous  form  of  our  surface  approximation  problem. 
We  formulated  it  in  the  required  abstract  form,  selected  appropriate  boundary  conditions,  and 
showed  that  it  is  well-posed  in  practice.  In  this  chapter,  we  face  the  task  of  applying  the  finite 
element  method  to  transform  the  variational  principle  into  an  appropriate  discrete  problem  whose 
discrete  solution  can  be  computed  fairly  easily.  Our  piecewise  continuous  representation  of  surfaces 
will  be  based  on  a  very  simple  finite  element  which  is,  howrver,  nonconforming.  This  will  force 
us  to  introduce  an  approximate  variational  principle  and  to  show  that  it  has  a  unique  solution 
which  converges  to  the  exact  solution  as  the  elements  shrink  in  size.  Only  then  can  we  undertake 
the  next  step  which  is  to  derive  the  discrete  problem  explicitly  as  a  linear  system  of  algebraic 
equations. 


18 


TERZOPOULOS 


MULTI-LEVEL  SURFACE  RECONSTRUCTION 


4.1.  Conforming  vs  Nonconforming  Methods 

Our  well-posed  variational  principle  satisfies  all  the  necessary  conditions  to  guarantee  that 
any  conforming  finite  element  method  applied  to  it  will  converge.  In  principle,  it  is  straightforward 
to  apply  a  conforming  finite  element  method  according  to  the  steps  in  Appendix  A.  We  generate 
a  finite  element  space  Sh  which  is  a  aubspace  of  our  admissible  space  V,  and  apply  the  Rits 
method  to  find  that  function  uh  6  SH  which  optimally  approximates  the  exact  solution  u  6  V. 
The  approximation  is  optimal  in  the  sense  that  it  is  closest  to  u  with  respect  to  the  strain 
energy  norm  a(-,  )^,  or  equivalently,  that  the  strain  energy  in  the  error  a(u  —  ufc,  u  —  tih)  is 
minimal.  To  construct  a  conforming  finite  element  subspace,  we  must  satisfy  the  completeness  and 
conformity  conditions  given  in  Section  A.4.  Since  the  energy  inner  product  a(-,  •)  contains  partial 
derivatives  of  order  m  —  2,  the  completeness  condition  requires  that  the  local  polynomial  defined 
within  each  element  subdomain  E  must  be  at  least  a  quadratic;  PE  D  n2(E)  for  all  E  6  TK 
(nn(f?)  denotes  the  set  of  ntb  degree  polynomials  over  E).  On  the  other  hand,  the  conformity 
condition  states  that  the  polynomials  must  be  of  class  Cl  across  inter-element  boundaries,  and 
consequently  Sh  C  C'(fi)  globally.  In  satisfying  both  conditions,  we  are  guaranteed  that  the  finite 
element  space  is  a  subspace  of  the  admissible  space  V,  and  that  there  exists  a  unique  optimal 
approximation  u k  £  Sk. 

If  fi  is  a  polygonal  region,  elements  with  straight  sides  will  suffice.  A  number  of  such  elements 
which  are  conforming  for  m  =  2  (i.e.,  problems  characterized  by  fourth-order  Euler-Lagrange 
equations)  are  available.  Examples  are  the  Argyris  triangle,  Bell  triangle,  and  Bogner-Foz- Schmidt 
rectangle  (see,  e.g.,  [Ciarlet,  1978],  [Strang  and  Fix,  1973],  [Zienkiewics,  1977]  and  the 
references  therein).  Unfortunately,  we  can  expect  serious  computational  difficulties  to  arise  in  the 
implementation  of  these  conforming  methods.  The  basic  source  of  difficulty  is  the  requirement 
of  continuity  of  first  partial  derivatives  across  inter-element  boundaries  —  either  the  structure 
of  the  conforming  element  spaces  PB  becomes  complicated,  or  their  dimension  is  large.  For  our 
problem,  the  simplest  conforming  polynomial  element  is  the  Bell  triangle,  in  which  we  have  a 
quintic  polynomial  uniquely  determined  by  18  nodal  variables  consisting  of  the  approximation  vk, 
as  well  as  its  first  and  second  partial  derivatives  at  the  three  vertices. 

As  is  described  in  Appendix  A,  the  dimensions  of  the  finite  element  space  can  be  reduced  by 
the  use  of  nonconforming  elements.  A  popular  nonconforming  element  for  fourth-order  problems 
is  Adini’t  rectangle,  whose  local  function  pE  is  a  12  degree-of-freedom  polynomial  with  nodal 
variables  being  the  approximating  function,  as  well  as  its  first  partial  derivatives  at  the  four 
vertices.  The  element  is  nonconforming  since  it  is  only  of  class  C°  across  inter-element  boundaries. 
Many  other  nonconforming  elements  have  been  developed  for  fourth-order  problems  (see,  e.g., 
[Ciarlet,  1978],  [Strang  and  Fix,  1973],  [Zienkiewics,  1977]). 

For  this  initial  implementation,  we  have  chosen  to  reduce  the  dimensions  of  the  finite  element 
space  as  much  as  possible  by  defining  what  for  our  problem  is  probably  the  simplest  successful 
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nonconforming  element  imaginable.  This  element  will  be  defined  next. 


4.2.  A  Simple  Nonconforming  Element 


We  will  define  a  finite  element  space  by  the  standard  procedure  outlined  in  Section  A.4. 
Suppose  that  ft  is  rectangular,  and  consider  a  uniform  triangutation  Th  of  Cl  into  identical  square 
elements  E,  where  the  fundamental  length  h  is  the  length  of  a  side  of  E.  By  definition,  we  require 
that  (Jeer*  E  —  Ci  and  that  the  elements  be  adjacent  and  overlap  along  their  sides.  A  point  in 
ft  is  a  node  of  the  triangulation  if  it  is  a  vertex  of  an  elemental  square  and,  as  usual,  we  consider 
the  elements  to  be  inter-connected  at  the  nodes.  The  nodal  variables,  will  simply  be  the  node 
displacements;  i.e.,  the  values  of  the  function  vh  £  Sh  at  the  nodes. 

The  next  step  is  to  define  a  space  PE  of  polynomials  pE  over  the  element  domain.  The 
polynomials  must  satisfy  the  completeness  condition  which  states  that  IT2  C  PE ,  since  the  energy 
inner  product  contains  derivatives  of  order  m  =  2.  This  is  the  requirement  that  the  polynomials 
be  able  to  reproduce  exactly  all  states  of  constant  strain  which,  in  this  case,  are  all  polynomials 
up  to  degree  two.  We  will  satisfy  this  requirement  by  choosing  PE  to  be  the  six-dimensional 
space  of  full  second  degree  polynomials  pE:  E  S?  such  that 


pE(x,  y)  =  ax 2  +  by2  -f  exy  +  dx  +  ey  +  /,  (1) 

where  the  six  real  parameters  a  bo  f  are  to  be  determined. 

We  must  ensure  that  pE  is  uniquely  specified  within  E  in  terms  of  the  node  displacements. 
To  do  so,  we  isolate  a  representative  element  and  set  the  origin  of  the  X-Y  coordinate  system 
at  its  lower  left  hand  corner,  as  illustrated  in  Figure  4.  Our  task  is  to  choose  a  pE -unisolvent 
set  of  nodes,  the  displacements  at  which  uniquely  determine  pE.  An  appropriate  choice  is 
the  six  nodes  shown  in  the  figure,  whose  node  displacements  are  denoted  by  v,tJ  £  SR,  for 
(i,j)  €  {( — 1,0),  (0, 0),  (1,0),  (0,  — 1),  (0, 1),  (1, 1)}.  Expressing  the  six  unknown  parameters  in  terms 
of  the  node  displacements  is  then  a  simple  matter  of  substituting  the  displacements  into  (1)  and 
solving  the  resulting  nonsingular  system  of  six  equations.  We  obtain 


°  =  2h2^Tl’°  —  2vo,o  "I"  v— i.o)> 

*  =  _  v°'- *)' 

c  =  '£a(vM  ~  vo, r  —  vi,o  +  v0,o), 

e=  2^(yo,i  -vo.-,), 

/  =  v0,o. 


(2) 


Of  course,  the  six  degrees  of  freedom  of  this  element  are  insufficient  to  enforce  C1  continuity 
of  vh  across  inter-element  boundaries.  Therefore,  the  element  is  nonconforming;  5*  £  C1  (ft).  It 
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Figure  4.  Uni»olvcnt  node*  for  the  nonconforming  clement. 


Y 


is  a  simple  matter  to  show  that  the  polynomials  pE  are  in  general  discontinuous  across  element 
boundaries,  although  continuity  is  maintained  at  the  nodes  themselves  because  each  polynomial 
interpolates  its  unisolvent  set  of  nodal  displacements.  At  this  point,  we  acknowledge  that  our 
element  is  somewhat  unorthodox  in  that  the  definition  of  pE  requires  nodal  variables  associated 
with  two  nodes  which  lie  outside  the  element  domain  E.  The  justification  for  this  transgression 
is  that  our  element,  as  defined,  will  yield  a  discrete  system  whose  matrix  is  particularly  simple 
and  uniform  in  structure.  This  will  simplify  the  eventual  implementation  considerably.  On  the 
other  hand,  alternate  arrangements  for  the  unisolvent  set  of  nodes  are  clearly  possible.  Perhaps 
a  more  appropriate  choice  from  a  biological  standpoint  would  be  a  hexagonal  .triangulation  with 
the  unisolvent  set  of  nodes  placed  at  the  vertices  of  hexagonal  element  domains  E  having  sides 
of  common  length  h.  Regardless  of  our  particular  choice,  the  quadratic  elements  must  first  be 
shown  to  be  convergent  since  they  will  invariably  be  nonconforming.  This  we  will  do  in  the  next 
section  for  the  simple  square  element. 
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4.3.  The  Approximate  Variational  Principle 


Due  to  the  nonconformity  of  the  elements,  Sh  <£_  l/2(fl)  =  V,  and  the  finite  element  space 
is  not  a  subspace  of  the  admissible  space  V.  Therefore,  in  lieu  of  the  energy  inner  product  a(-,  ■) 
of  equation  (3.3),  we  are  forced  to  substitute  the  approximate  energy  inner  product  a/,(-,  •),  given 
by  (A. IS),  which  can  be  written  as 


»h)  =  E  /  /  AtsfcAvfc  —  (1  —  *)(«£,  +  is^w^  -2uhIvvKzy)dxdy 

£eT*  E 

+  0  E  u/l(*..y.)vh(*..v.) 

(i..v.)ec 

=  E  /  L  A“fcA«*  +  (1  -  +  2uhzvv*v  +  «$„«£„)  dx  dy 

eer»  E 


+  0  E  wh(*«.V.)wh(*..V«). 

(*.,»<  )€C 


(3) 


where  0  <  a  <  1.  The  corresponding  approximate  variational  principle  and  variational  equation 
are  given  by  equations  (A. 17)  and  (A. 18)  respectively. 

Does  the  approximate  variational  principle  have  a  unique  solution  uh  6  S*1?  To  answer  this 
question,  we  proceed  in  the  spirit  of  Section  A. 5  by  equipping  Sh  with  a  norm  which  we  will 
employ  to  show  that  ah(-,  ■)  is  uniformly  S '‘-elliptic. 

Proposition  3.  The  mapping  ||v'l||/l:  Sk  >-»  5R  defined  by 

E  l*fcll.K+  E 

v£eT<-  (x„Vi  )ec 


where  \vh\ 2,e  =  (/  Ie(vzx)2  +  (uiy)2  "+•  (vyV)2  dxdy)^  is  the  second-order  Sobolev  semi-norm  (see 
(A. 3)),  is  a  norm  over  Sh. 

Proof,  n  il  is  a  priori  only  a  semi-norm  over  Sh.  Consider  a  vh  £  Sh  such  that  ||t>h||fc  =  0. 
Then  it  must  be  the  case  that  (i)  |vh|2,£:  =  0  for  ail  E  £  Th,  and  that  (ii)  vh(x,,y,)  =  0  for 
all  (x,,i/i)€  C.  Because  of  their  interpolatory  nature,  the  local  polynomials  pE  are  continuous  at 
all  the  nodes.  Moreover,  by  condition  (i),  vh  must  be  a  first-degree  polynomial  inside  every  E. 
With  a  =  6  =  c=  0in  equation  (2),  it  is  a  simple  matteT  to  show  that  this  implies  that  vh  is  a 
continuous  linear  function  over  fi.  Now,  by  condition  (ii),  vh  is  sero  at  all  (x,,  y, )  £  C.  Since  the 
continuous  problem  is  assumed  to  be  well-posed,  C  contains  at  least  three  noncoilinear  points. 
Consequently  vh  =  0,  and  ||-||k  is  therefore  a  norm.  | 

Proposition  4.  The  approximate  energy  inner  product  ah(-,  •)  is  uniformly  S '‘-elliptic. 

Proof. 
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<l/.(v\  t/“)  =  Yl  //.a(Awh)2  -t  t1  “^(K'x)2  %»xy)2  +{VyvY]<tody  +  0  E 

EeT’-  (x..v.)ec 

=  0-*)(  E  f  [(W +  (<)*  + K)2d**y  +  E  «*(*..  v.)2) 

yiceT’>J  Jk  (i„v.)ec  ' 

+  E  f  [  <7(Avh)2  +  (1  -  ff)(viy)2  dxdy +  (^  +  cr  —  1)  ]T  w'‘(ilt»,)2 

t  GT'-'  ■/k 

>  (1  -<0||v'1||£,  for  0<ff<l,  P  >  l  —  a. 

Since  1  — -  cr  is  positive  for  0  <  a  <  1,  o^(-,  •)  is  uniformly  Sfc-elliptic. 


K.Vilec 


Therefore,  the  approximate  variational  principle  has  a  unique  solution  uh  £  Sh .  Moreover, 
because  the  ellipticity  is  uniform,  Strang’s  lemma  (Theorem  A. 5)  applies,  and  uh  will  converge 
to  the  exact  solution  u  6  V  as  h  — *  0,  if  the  approximation  is  consistent  in  the  sense  of  equation 
(A. 21).  To  verify’ consistency,  we  apply  the  patch  test,  Theorem  A. 6. 


Proposition  5.  The  square,  nonconforming  element  whose  local,  quadratic  function  iB  defined  by 
(1)  and  (2)  passes  the  patch  test. 

Proof.  Consider  an  arbitrary  patch  of  four  adjacent  elements,  all  of  which  share  a  common  node 
v,^  internal  to  the  patch,  as  shown  in  Figure  5.  Now,  suppose  that  we  impose  a  constant  strain 
condition  on  the  patch;  that  is,  suppose  that  we  constrain  the  displacements  at  all  remaining 
nodes  around  the  periphery  of  the  pr  tch  by  assigning  to  them  values  consistent  with  the  function 
ir2  £  I!2,  an  arbitrary  second-degree  polynomial.  Next,  we  solve  the  approximate  variational 
principle  (A. 17)  over  the  patch  domain.  This  reduces  to  solving  for  the  unknown  displacement 
at  the  common  unconstrained  node  vtiJ  such  that  it  minimises  a  quadratic  equation.  It  is  a 
matter  of  routine  algebraic  manipulation  to  show  that  the  displacement  obtained  will  also  be 
consistent  with  ir2  (we  omit  the  details).  In  fact,  one  can  show  that  this  is  true  for  the  internal, 
unconstrained  nodes  of  an  arbitrary  patch  of  any  number  of  elements  whose  boundary  nodes  are 
made  consistent  with  7t2.  Therefore,  the  element  passes  the  patch  test.  | 


Having  proved  the  above  propositions,  we  can  now  be  secure  in  the  fact  that  our  approximate 
variational  principle  will  provide  unique  discrete  solutions  which  will  converge  to  the  exact  solution 
of  the  continuous  problem  as  the  discretization  is  made  increasingly  fine.  A  limit  to  the  order  of 
convergence  that  we  can  expect  from  our  approximation  is  given  by  (A. 15)  —  since  our  element 
is  complete  only  through  quadratics  (fc  =  2),  we  are  limited  to  a  convergence  rate  of  order  h2  in 
displacement  (a  =  0).  For  a  more  precise  statement,  we  should  take  into  account  the  consistency 
error  term  in  equation  (A. 20).  Nevertheless,  we  will  bypass  this  complicated  analysis  because  the 
consistency  error  is  not  expected  to  be  large  for  smooth  u  which  is  normally  the  case  when 
approximating  smooth  surfaces. 

4.4.  The  Discrete  Problem 

We  are  finally  ready  to  derive  an  explicit  form  of  the  discrete  problem  associated  with  our 
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Figure  5.  Applying  the  patch  test  to  four  adjacent  element*. 


i  1 


approximate  variational  principle.  There  are  essentially  two  ways  of  proceeding.  One  possibility  is 
to  find  the  Ritz  basis  functions  <f>,  which  are  associated  with  our  finite  clement  and  which  span 
the  space  Sk.  The  basis  functions  are  nonconforming  piecewise  quadratics  with  local  support,  and 
a  basis  function  associated  with  each  node  of  the  triangulation.  We  can  then  use  the  variational 
equation  (A. 18)  directly,  and  write  the  discrete  problem  as  the  linear  system  of  equations  analogous 
to  equation  (A.13)  by  computing  the  matrix  coefficients  4>.).  Unfortunately,  the  piecewise 

continuous  "nature  of  the  basis  functions  makes  them  tedious  to  manipulate,  especially  near 
the  boundary.  We  will  adopt  an  alternate  approach  which  altogether  avoids  the  derivation  and 
manipulation  of  the  basis  functions.  The  approach  is  to  solve  the  approximate  variational  principle 
by  minimizing  the  functional  £),(•)  of  equation  (A. 17).  Before  doing  so,  however,  we  make  two 
additional  simplifications. 

The  first  simplification  involves  taking  a  conservative  stance  once  more.  There  is  no  reason 
to  believe  that  the  human  visual  system  is  biased  in  the  depth  values  it  assigns  such  as,  for 
example,  making  all  of  them  too  small  or  too  large.  That  is  to  say,  we  have  no  reason  to  assume 
that  there  is  an  external  influence  on  the  surface  other  than  that  provided  by  the  constraints  C, 
and  we  should  therefore  nullify  the  externally-applied  surface  force:  g[x,  y)  =  0.  The  linear  form 
(3.4)  for  the  free  plate  then  reduces  to 

f(vh)=0  Y,  (cl*»y.>wh(*«.V«)  —  JC(,€,*)Y  (4) 

(n,v.)€C  '  ' 
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The  second  simplification  involves  the  choice  of  a  numerical  value  for  the  constant  a  in  our 
approximate  energy  inner  product  ah(-,  •)  given  by  (3).  According  to  the  proof  of  Proposition  4, 
we  are  at  liberty  to  choose  any  value  in  the  range  0  <  a  <  1,  therefore,  the  simplifying  choice 
<7  =  0  will  be  made.4  Setting  <7  =  0  in  (3),  the  energy  inner  product  becomes 

«*(*\**)  =  £  f  fyiyiI+2uhiyzv+^yvydxdy  +  fi  £  uh(zt,y,)vh(z„y).  (5) 
E€ T k  J  JE  (xovoec 

Thus,  according  to  equation  (A. 17),  we  obtain  the  simplified  energy  functional 

^ah(v\vE)~f(vh) 

=  5  £  /  /(*4)2  +  2(*4 )2  +  K„)2didy +  |  £  [vh(x,,  y,)  c(x,,y,)]‘ •  ^ 

z  e& r*J  JE  ^(x<,y<)ec 

The  expression  inside  the  integral  will  be  recognized  as  the  " quadratic  variation"  expression  used 
by  Grimson  [1981a]. 

Since  the  triangulation  over  the  rectangular  region  fi  is  that  of  a  uniform  square  grid,  it  is 
convenient  to  impose  on  the  nodes  the  natural  lexicographic  indexing  scheme  implied  by  Figure 
4.  We  index  the  nodes  by  («,j)  for  i  =  1  ,...,NX  and  j  =  \,...,NV  where  N7  and  Nv  are  the 
number  of  nodeB  along  the  x  and  y  axis  respectively.  The  total  number  of  nodes  is  N  =  Nz  X  Ny. 
The  displacement  at  node  (i,j)  is  denoted  by  the  variable  vj1^,  and  all  the  displacements  together 
are  denoted  by  the  vector  vh  £  &N . 

The  next  step  is  to  express  the  functional  in  terms  of  the  node  displacements  with  the  help 
of  our  element.  Inside  each  element  domain  E,  uh  is  a  quadratic  polynomial  given  by  (1)  and  (2). 
Therefore,  the  second  partial  derivatives  of  pE  are  con.-1’  -  wit  within  E,  and  are  given  by 

*4 Is  =  pfx  =  2o  =  -  2ve}  +  v?_it>); 

vyV\ E  =  Pyy  =  26  =  +  I  “  2v<4  +▼<*.,- 1)> 

vxy\ £  =  Pfy  =  c  =  ^(▼h-i.,  +  1 

where  it  is  assumed  that  i,j  is  the  index  of  the  lower  left  hand  node  of  E.  The  form  of  these 
second  derivatives  will  be  recognized  as  being  simply  the  finite  difference  approximations  of  order 
h 2  for  the  respective  derivatives  on  a  uniform,  square  mesh  (see,  e.g.,  [Abramowit*  and  Stegun, 
1965,  pg.  884]).  Of  course,  this  result  is  a  consequence  of  our  particular  choice  of  finite  element, 
and  it  will  lead  to  a  particularly  simple  discrete  problem.  With  other  elements  one  cannot  expect 
to  obtain  finite  difference  expressions,  even  for  uniform  triangulations.  Substituting  the  expressions 

4Recall  that  o  is  the  Poisson  constant  of  the  elastic  material,  so  our  choice  implies  that  the  material 
does  not  change  in  width  at  it  stretches  lengthwise.  Although  this  value  is  not  meaningful  physically,  it 
is  perfectly  acceptable  mathematically.  Aside  from  a  question  of  convenience,  there  is  further  evidence 
that  supports  this  choice  in  terms  of  the  optical  laws  of  image  formation  (see  [Grimson,  1982b]). 
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for  the  derivative!  into  a^v*1,  vh)  in  (6)  and  noting  that  the  area  of  each  element  is  h2,  we 
obtain5 


We  can  write  the  above  expression  for  the  functional  (aside  from  the  additive  constant  term) 
in  matrix  form  as 

ffc(vk)-  i(v'*,AV)-(f\v'‘),  (7) 

where  (•,  ):!RN  X  >-»  !R  denotes  the  familiar  Euclidean  inner  product,  and  Ah  £  iHNN  is  a 
matrix  of  coefficients.  Clearly,  equation  (7)  is  the  discrete  equivalent  of  the  functional  (6).  For 
the  linear  term,  we  have  f*  =  0ch  where  eh  £  is  a  vector  whose  entries  associated  with 
constrained  displacements  are  the  constraint  values  c* and  the  rest  are  zero.  On  the  other  hand, 
the  matrix  Ah  £  5RNN  which  forms  the  quadratic  term,  is  a  matrix  of  coefficients  which  can 
be  brok  :n  down  as  the  sum  of  two  matrices:  Ah  =  Aj,  -(-  Bh.  The  matrix  Bh  is  a  diagonal 
matrix  whose  diagonal  entries  associated  with  constrained  displacements  are  equal  to  0 ,  and  the 
remainder  are  sero.  As  is  clear  from  equation  (A. 13),  the  entries  of  the  other  component  Aj  can 
be  interpreted  as  inner  products  between  pairs  of  basis  functions  of  the  finite  element  space  Sh. 
Since  the  basis  functions  have  local  support,  most  of  these  inner  products  will  be  zero  thereby 
making  Ah  sparse  and  banded.  Moreover,  since  by  Propositions  1  and  4,  the  energy  inner  product 
is  symmetric  and  S^-elliptic,  Ah  is  a  positive  definite,  symmetric  matrix.  These  are  important 
properties  from  a  computational  point  of  view. 

To  obtain  the  minimum  of  Eh('/h)  we  set  to  zero  its  partial  derivatives  with  respect  to  each 
of  the  displacements  .  The  minimizing  vector  of  displacements  uh  satisfies  the  condition 

v£*(u'1)  =  aV1  -  r*  =  #, 

(where  V  is  the  gradient  operator)  which  yields  a  discrete  problem  in  the  form  of  a  system  of 
linear  equations: 

A'V  =  f*,  (8) 

where  the  entries  of  Ah  are  given  by 

5We  also  assume  for  simplicity  that  all  constraints  coincide  with  nodes  (t, j)  in  Th.  Hence,  we 

will  denote  the  constraints  by  c* 
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omit  the  details  and  give  the  final  result  in  terms  of  a  set  of  computational  molecules  which  are 
illustrated  in  Figure  6  in  relation  to  the  lower  left  hand  edge  of  Cl  whose  boundary  is  indicated  by 
bold  links.  Obviously,  computational  molecules  for  the  remaining  edges  are  appropriate  rotations 
of  those  shown.  The  particular  computational  molecule  associated  with  a  node  specifies  the 
(nonzero)  coefficients  of  the  equation  for  that  node.  For  example,  the  equation  for  the  displacement 
at  a  node  (i,j)  in  the  interior  of  the  region  (indicated  by  the  double  circle  in  the  topmost 
computational  molecule  in  Figure  6)  is 

-  tfW- W  +  +  <>-1  +  <r+i) 

+  nW-M-l  +  1  +  +  vH-i.j  +  i) 

1  (9) 

+  +  VH-2,>  +  Vi.,-2  +  ^.j+a) 

+  5^  +  K,  - 

The  terms  involving  (3  are  present  only  if  there  is  a  constraint  c,,;  at  that  node. 

The  sparseness  of  \h  is  evident  from  the  above  equation  —  matrix  rows  associated  with 
interior  nodes  have  only  13  nonzero  entries,  while  rows  associated  with  nodes  near  the  boundaries 
have  even  fewer.  Also,  note  that  the  computational  molecule  for  the  center  of  the  region  is  a 
factor  ct  h'2  (due  to  the  elemental  area)  times  the  Unite  difference  approximation  of  order  h2 
for  the  biharmonic  operator  [Abramowilz  and  Slegun,  1965,  pg.  885]  that  is  associated  with  the 
Euler- Lagrange  equation  for  our  variational  principle  This  is  an  expected  consequence  of  our 
particular  choice  of  element  which  yielded  finite  difference  app-oximations  for  the  second  partial 
derivatives  of  v1’.  Moreover,  aside  from  multiplicative  constants,  the  same  molecules  were  obtained 
by  Grimson  [1981a]  in  the  specification  of  a  (conjugate  gradient)  mathematical  programming 
algorithm.  As  was  previously  argued  however,  the  fini*c  element  method  is  richer  in  that  it 
systematically  suggests  many  alternative,  less- restrictive  inangulations,  as  well  as  more  general 
local  representations  for  surfaces. 

5.  MULTI-LEVEL  SURFACE  RECONSTRUCTION 

As  we  have  seen,  the  application  of  the  finite  element  method  to  a  well-posed  quadratic 
variational  principle,  such  as  the  one  on  which  our  surface  approximation  problem  is  based, 
inevitably  leads  to  an  equivalent  discrete  problem  which  takes  the  form  of  a  linear  system  of 
algebraic  equations.  The  matrix  of  coefficients  of  this  nonsingular  system  is  symmetric,  positive 
definite,  sparse,  and  banded.  Computing  the  most  consistent  approximating  surface  now  amounts 
to  solving  this  system  and,  in  this  chapter,  we  adopt  an  efficient  hierarchical  algorithm  to  perform 
this  task.  We  will  proceed  to  develop  the  algorithm  and  to  demonstrate  its  performance  it  on 
constraints  from  stereopsis. 
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5.1.  Possible  Methodologies  for  Solving  the  Discrete  Problem 

The  solution  of  linear  systems  is  a  very  important  problem  in  numerical  analysis  and  the 
many  techniques  which  have  been  developed  fall  into  essentially  two  broad  classes  —  direct 
methods  which  yield  the  solution  in  a  finite  number  of  steps,  and  iterative  methods  which  typically 
converge  asymptotically  to  the  solution  (see,  e.g.  [Dahlquist  and  Bjorck,  1974]  or  [Glad well  and 
Wait,  1979]). 

Direct  methods  include  matrix  inversion  methods  such  as  Gaussian  elimination  and  LU 
decomposition.  Although  widely  used  for  solving  finite  element  equations,  they  usually  do  rot 
exploit  the  sparseness  and  bandedness  of  the  system  matrix  because,  during  the  inversion  process, 
the  sparse  matrix  is  transformed  into  a  full  one.6  Consequently,  all  the  elements  of  the  matrix 
must  be  stored.  Moreover,  direct  methods  are  typically  global  and  sequential  algorithms,  which 
makes  them  unsuitable  models  for  neurally-based  visual  processes. 

On  the  other  hand,  the  class  of  iterative  methods  readily  gives  rise  to  biologically-fe&sible 
algorithms.  Examples  in  this  class  are  relaxation  methods  such  as  Jacobi  relaxation,  Gauss-Seidet 
relaxation  and  successive  overrelaxation,  as  well  as  gradient  methods  such  as  gradient  descent  and 
the  conjugate  gradient  method.  Iterative  methods  exploit  the  sparseness  of  the  matrix  inasmuch 
as  they  do  not  modify  its  elements  from  one  iteration  to  the  next.  Therefore,  only  the  relatively 
few  nonzero  matrix  elements  need  be  stored.  Owing  to  the  sparseness  and  banded  structure 
of  the  matrix,  iterative  methods  require  local-support  computations  only,  and  in  certain  forms 
such  as  Jacobi  relaxation  and  gradient  methods  the  computations  can  be  performed  in  parallel. 
Because  iterative  methods  in  general  and  relaxation  in  particular  are  fundamental  to  the  ensuing 
discussion,  an  introduction  to  some  of  the  relevant  mathematics  of  this  class  of  techniques  is 
included  in  Appendix  B  for  convenience. 

The  algorithms  we  are  contemplating  are  to  be  executed  by  computational  mechanisms  in 
the  form  of  networks  of  many  simple  processors,  such  as  neurons,  which  are  directly  connected 
only  to  near  neighbors.  Due  to  the  myopic  nature  of  the  processors,  global  interactions  can  take 
place  only  indirectly,  through  the  iterative  process,  by  an  incremental  propagation  of  information. 
Normally,  the  network  is  large  and  since  this  is  reflected  in  the  size  of  the  linear  system,  we 
anticipate  that  a  vast  number  of  iterations  will  be  required  for  any  relaxation  or  gradient  method 
to  converge.  Typically,  the  number  of  iterations  will  be  on  the  order  of  Nm,  where  N  is  the 
dimension  of  the  matrix,  and  m  is  the  highest  order  of  partial  derivatives  present  in  the  energy 
inner  product,  which  in  our  case  is  two.  Grimson’s  [1981a]  formulation  of  surface  interpolation  as 
a  problem  in  mathematical  programming  naturally  led  him  to  the  choice  of  a  gradient  method 
for  its  solution  and,  not  unexpectedly,  disappointingly  slow  convergence  rates  were  observed  due 
to  the  large  size  of  the  images  typically  encountered. 

*For  a  positive  definite  symmetric  matrix,  the  inverse  matrix  remains  banded,  but  is  no  longer  sparse 
within  the  band.  The  inverse  matrix  is  the  discrete  Green’s  function  for  our  problem,  which  in  general 
has  global  support  over  ft. 
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Recently,  a  class  of  iterative  techniques  called  multi-grid  methods  have  seen  increased 
application  to  the  numerical  solution  of  boundary  value  problems  for  which  they  achieve 
convergence  in  essentially  order  N  number  of  operations  (Brandt,  1977a,  1977b;  Brandt  and  Dinar, 
1979].  This  spectacular  improvement  results  from  the  use  of  a  hierarchy  of  grids  to  increase  the 
efficiency  of  the  global  propagation  of  information.  Multi-grid  algorithms  feature  both  intra-grid 
and  inter-grid  processes.  Typically,  the  intra-grid  processes  are  relaxation  iterations,  while  the 
inter-grid  processes  are  local  polynomial  interpolations.  Therefore  the  multi- grid  algorithms  are, 
in  principle,  biologically  feasible.  A  final  issue  which  speaks  in  favor  of  adopting  them  to  vision 
is  the  intrinsic  multi-level  structure  of  the  earliest  stages  of  the  visual  system  itself  and,  as  we 
argued  in  the  introduction,  the  apparent  need  to  maintain  this  structure  at  least  to  the  level  of 
the  2  J-D  sketch. 

We  therefore  advocate  a  hierarchical  approach  to  surface  reconstruction,  which  we  will 
develop  initially  in  the  context  of  the  Marr-Poggio  stereo  theory  whose  clear  multi-level  structure 
provides  ample  motivation.  At  the  heart  of  the  proposed  scheme  lies  a  multi-grid  algorithm 
adapted  to  the  fast  solution  of  a  hierarchy  of  discrete  thin  plate  surface  approximation  problems. 
In  the  following  sections,  we  present  the  underlying  theory  and  build  up  a  detailed  description  of 
the  multi-level  algorithm. 

5.2.  The  Multi-Level  Equations 

As  we  have  stated,  the  stereo  module  generates  sparse  depth  information  over  a  range  of 
resolutions.  The  information  at  any  particular  scale  can  be  thought  of  as  a  set  of  constraints 
which,  at  that  level,  define  a  well-posed,  discrete  surface  approximation  problem.  It  is  natural 
then  to  view  our  surface  reconstruction  problem  as  the  solution  of  a  hierarchy  of  such  discrete 
problems.  The  discretizations  axe  performed  in  the  usual  way  by  introducing  a  sequence  of  finite 
element  spaces  Sh' , . . . ,  Shl-  over  the  rectangular  domain  ft,  where  L  is  the  number  of  levels 
and  hi  >  •  •  •  >  ht  are  the  fundamental  lengths  of  the  elements  at  each  level.  In  the  familiar 
notation,  we  will  denote  the  functions  which  are  members  of  the  finite  element  spaces  by  (italic 
face)  vllk  6  Shk,  and  the  parameters  (i.e.,  the  nodal  displacements)  which  define  these  functions 
according  to  (A. 11)  by  (bold  face)  vectors  vhk  6  ,  where  Nhk  is  the  dimension  of  Shk .  The 

hierarchy  of  problems  is  then  given  by  the  sequence  of  L  linear  systems  (see  equation  (4.8))  of 
the  form 

A'“u'“  1  <k<L,  (1) 

whose  discrete  solutions  uh*  €  for  1  <  k  <  L  define  a  sequence  of  functions  uhk  6  Sk  which 

constitute  the  hierarchy  of  full  surface  representations. 

Although,  in  theory,  there  need  be  no  restriction  in  the  relationship  of  element  sizes  from 
one  level  to  the  next,  a  number  of  practical  implementation-related  issues  point  towards  the 
subdivision  of  each  square  element  domain  on  a  given  level  Sh>  into  four  identical  element  domains 
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on  the  next  finer  level  Shk  +  l;  that  is,  we  choose  h *  =  2h/,+i.  Consequently,  Shk  wi!l  be  a 
subspace  of  S,l*+1,  and  the  implementation  of  the  inter-level  processes  is  simplified  substantially. 
Moreover,  the  2:1  ratio  is  a  natural  one  in  view  of  the  spatial- frequency  bandpass  channels  in  early 
vision  whose  center  frequencies  are  spaced  approximately  onts  octave  apart,  the  spatial  resolution 
of  a  channel  being  about  twice  that  of  the  immediately  coarser  one  [Wilson  and  Giese,  1979]. 
Finally,  the  choice  can  be  shown  to  be  near  optimal  in  terms  of  the  multi-grid  convergence  rate 
[Brandt,  1977a,  pg.  353].  Since  the  triangulation  of  Cl  associated  with  our  simple  elements  is  a 
uniform  grid  of  square  element  domains,  the  2:1  ratio  implies  that  in  scanning  along  the  x  or 
y  directions,  every  second  node  of  a  grid  coincides  with  a  node  of  the  next  coarser  grid  and, 
furthermore,  that  the  number  of  nodes  is  related  from  one  level  to  the  next  by  Nhk-'  —  $Nhi. 
Therefore,  the  total  amount  of  space  required  to  maintain  all  of  the  representations  is  bounded 
by  JVhL(  1  +  i  4;  +  •  ■  •)  =  i N1*1--,  i.e.,  it  is  only  a  small  fraction  more  than  that  required  for 

the  finest  grid. 

One  can  think  of  several  possibilities  for  exploiting  the  hierarchy  of  discrete  problems  to 
increase  the  convergence  rate  of  the  iterative  process.  Perhaps  the  first  idea  that  comes  to  mind 
is  to  solve  the  system  at  the  coarsest  level,  which  can  be  done  very  quickly,  and  use  the  solution 
as  an  initial  approximation  in  the  iterative  solution  of  the  next  finer  level,  proceeding  in  this 
manner  to  the  finest  level.  This  is  an  effective  acceleration  strategy  which  is  almost  as  old  as  the 
idea  of  relaxation  itself  [Southwell,  1946].  Although  it  is  suitable  for  obtaining  a  single  accurate 
solution  at  the  finest  level,  it  cannot  generate  solutions  having  the  finest-level  accuracy  over  the 
hierarchy  of  coarser  levels,  since  the  approximation  error  increases  as  the  elements  become  larger. 
This  is  undesirable  from  the  point  of  view  of  our  surface  reconstruction  problem.  Here  we  require 
that  the  accuracy  of  the  finest^resolution  surface  be  maintained  throughout  the  coarser  surface 
descriptions.  This  will  guarantee  that  the  shape  of  the  surface  will  be  consistent  over  the  hierarchy 
of  representations. 

The  stipulation  that  accuracy  be  maintained  is  further  motivated  by  psychophysical  studies 
into  the  phenomenon  of  visual  hyper  acuity  (see,  e.g.,  [Westheimer,  1977;  Westheimer  and  McKee, 
1975,  1977]).  Related  computational  studies  indicate  that,  in  principle,  sharp,  well-defined  intensity 
edges  can  be  localized  to  high  (sub-receptor  separation)  accuracies  from  the  V2G  convolution 
values  through  a  process  of  spatiotemporal  interpolation  [Barlow,  1979;  Marr  et  al.,  1980;  Hildreth, 
1980;  Crick  et  al.,  1981;  Fahle  and  Poggio,  1981].  Consequently,  it  seems  that  although  the  depth 
constraints  arising  from  the  larger  channels  in  stereopsis  represent  coarser  spatial  samplings  of  the 
scene,  excluding  erroneous  matches,  the  samples  may  provide  highly  accurate  range  information. 

The  only  way  that  consistent  accuracy  can  be  maintained  throughout  the  hierarchy  of  full 
surface  representations  is  to  allow  the  coarser  levels  access  to  the  high- resolution  information  in 
the  finer  levels.  The  multi-grid  algorithm  provides  such  a  flow.  The  hierareny  of  levels  cooperate, 
through  a  bi-directional  flow  of  information,  to  simultaneously  generate  multiple,  equally- accurate 
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surface  representations,  and  do  so  with  much  less  computational  effort  than  would  be  expended 
in  solving  the  finest-level  system  in  isolation.  To  see  how  this  is  accomplished,  we  will  initially 
consider  only  two  levels,  a  fine  one  and  a  coarse  one,  associated  with  the  finite  element  spaces  Sh 
and  S2h  respectively.  Suppose  that  by  some  iterative  process  we  obtain  an  approximate  solution 
v2h  to  the  coarse  level  system  A2/lu2h  =  fah,  which  is  then  interpolated7  to  the  fine  level  where 
it  becomes  the  initial  approximation  vh: 

yh  4—  I  y2\ 

2h=*/i 

The  mapping  hh=*h  '■  S2h  >-*  Sh  denotes  interpolation  from  the  coarse  space  to  the  fine  space. 
Normally,  \h  will  require  substantial  improvement. 

Let  uh  be  the  solution  to  the  fine  level  system  \huh  —  f*.  Then  we  can  define  the  error 
vector  in  a  given  approximation  vh  by  eK  —  u^  —  vh.  Clearly,  if  eh  could  be  computed,  it  could  be 
added  as  a  correction  to  vh,  thereby  giving  us  the  desired  solution.  But  because  the  computation 
of  e'1  would  take  about  as  much  effort  as  computing  uh  itself,  doing  so  would  not  be  helpful. 
On  the  other  hand,  if  we  could  somehow  approximate  the  error  function  fh  by  a  function  e2h  in 
the  coaise  space  S2h ,  such  an  approximation  can  be  obtained  quickly  due  to  the  fact,  that  the 
coarse  space  has  only  one  quarter  the  dimensionality  of  the  fine  space.  Such  an  approximation  is 
generally  not  possible,  however,  because  e'1,  having  been  generated  by  an  interpolation  from  the 
coarse  grid  solution,  is  certain  to  have  large  fluctuations  with  wavelengths  less  than  Ah.  These 
high-frequency  Fourier  components  could  not  be  approximated  on  the  coarse  grid  because  there 
they  would  alias  as  lower-frequency  components.  Before  a  meaningful  approximation  to  the  error 
can  be  obtained  on  the  coarse  grid,  the  high-frequency  components  must  be  eliminated;  that  is 
to  say,  the  error  function  eh  must  be  smoothed. 

Since  smoothing  is  inherently  a  local  operation,  it  should  not  be  surprising  that  local  iterative 
methods,  inefficient  as  they  are  in  obtaining  solutions,  are  very  efficient  at  smoothing  the  error 
function.  In  particular,  although  relaxation  generally  requires  very  many  iterations  to  eliminate 
the  global,  low-frequency  components  of  the  error,  it  only  takes  a  few  iterations  to  eliminate 
the  local,  high-frequency  components.  This  behavior  can  be  predicted  mathematically  by  a  local 
Fourier  analysts  of  the  given  iterative  method  [Brandt,  1977a].  The  analysis  involves  a  local  Fourier 
expansion  of  the  error  function  followed  by  an  examination  of  the  effect  that  a  single  iteration 
has  on  the  amplitudes  of  each  component.  An  important  quantity  which  is  obtained  through  this 
analysis  is  the  smoothing  factor  p  of  the  iterative  scheme,  which  is  defined  as  the  worst  (i.e., 
the  largest)  amplification  of  a  high-frequency  component  of  the  error  from  one  iteration  to  the 
next.  As  an  example,  in  Appendix  C  we  carry  out  a  local  Fourier  analysis  of  the  appropriate 
Gauss-Seidel  scheme  for  our  discrete  problem,  and  show  that  p  «  0.8.  This  implies  that,  for  our 
problem,  ten  Gauss-Seidel  iterations  on  the  fine  grid  are  sufficient  to  reduce  the  high-frequency 

7  Lagrange  interpolation  of  a  suitable  order  may  be  used. 
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components  of  eh  by  approximately  an  order  of  magnitude.  A  more  effective  weighted  Jacobi 
relaxation  scheme,  which  is  also  suitable  for  our  problem  and  gives  a  p  —  0.549,  is  described  in 
[Brandt,  1977a,  pg.  342].8 

Once  the  error  has  been  smoothed,  it  may  be  inexpensively  approximated  on  the  coarse 
grid.  The  equation  for  eh  on  the  fine  grid  is  the  residual  equation 

Aheh  =  r\  where  r'*  =  I*  -  Ahvh  (3) 

is  called  the  residual  of  the  approximation  vh.  The  approximation  to  this  equation  on  the  coarser 
grid  is 

A2ke2h  =  J 

h=*2A 

where  the  mapping  I/i=,2a  :  Sh  i->  SiK  is  an  “interpolation”  from  the  fine  space  to  the  coarse 
space.  Because  S2h  C  Sh,  the  mapping  can  be  a  simple  injection  or  some  local  averaging  of  node 
displacements  from  the  fine  grid  to  the  coarse.  After  e2h  is  computed,  a  better  approximation  to 
vh  may  be  obtained  by  interpolating  the  coarse  grid  correction  back  to  the  fine  grid;  that  is,  by 
making  the  replacement: 

\h  «-  yK  +  I  e2\ 

2/i=*h. 

This  correction  practically  annihilates  the  smooth  part  of  the  error  eh. 

Brandt  [1977a,  1977b]  calls  the  foregoing  scheme  a  correction  scheme  in  view  of  the  fact 
that  the  function  computed  on  the  coarse  grid  is  the  error  function;  that  is,  the  correction  to 
the  fine  grid  approximation.  The  correction  scheme  is  easy  to  implement,  but  it  is  unsuitable  for 
our  surface  reconstruction  problem  because  instead  of  an  error  function  e2h,  we  require  that  the 
function  computed  in  the  coarse  space  be  a  function  u2h  which  represents  explicitly  the  distances 
to  surfaces  in  the  scene.  This  may  be  accomplished  by  a  reformulation  of  the  correction  scheme 
equations  which  converts  them  into  those  of  the  related  full  approximation  scheme. 

First,  we  rewrite  (3)  in  the  equivalent  form 

AV  +  eh)-AV  =  A 

which  may  be  approximated  by  the  corresponding  coarse-grid  equation 

A2h{  j  v*+e2A)_A2h(  I  yh)  =  I  f\ 
h—2h  A-*2h  h-»2h 

Defining  a  new  function  u2h  in  Slh  by  the  nodal  displacement  vector  u,h  —  Ih-*aa  vh  -f-  e2h,  we 
obtain  the  coarse  level  system 

'Brandt  proposes  this  scheme  for  solving  biharmonic  boundary  value  problems.  The  scheme  is  also 
appropriate  for  our  surface  approximation  problem  which  is  in  essence  a  biharmonic  problem  in  view  of 
the  associated  Euler-Lagrange  equation. 
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AJ'*u2'*  =  g2H,  where  g: 


2fc 


=  \2h{  I  vh)  +  I  r* 
k-*2k  k=*2k 


(4) 


It  is  natural  to  interpret  (4)  as  the  original  coarse-level  system  A2hu2h  =  but  with  a  right-hand 
side  which  has  been  modified  using  information  from  the  fine  grid  so  as  to  maintain  the  fine-grid 
accuracy  in  the  coarse-grid  function  u2h .  Thus,  g"'h  is  an  estimate  of  the  local  truncation  error  on 
the  coarse  level  relative  to  the  fine  level  (see  [Brandt,  1977b,  pg.  284]). 

Once  the  solution  u2h  of  equation  (4)  is  available,  we  can  write  e2h  =  u2'“  —  Ik=»2k  vh  as  the 
desired  coarse-grid  approximation  to  the  fine-grid  error,  and  the  approximation  on  the  fine  level 
can  be  corrected  by  the  replacement 


vh  «-  v*  +  I  (u2*  -  I  v*). 

2k-»k  k=-»2k 


(5) 


Note  that  since  l2k=*k  Ik=»2h  vh  vh  the  replacements  given  by  (2)  and  (5)  are  not  equivalent. 
Since  u2A  is  a  low- frequency  correction,  the  replacement  indicated  by  (2)  would  destroy  the 
high-frequency  components  of  vh  whereas  the  replacement  indicated  by  (5)  preserves  them. 


How  do  we  solve  the  coarse-grid  equation  (4)?  The  obvious  answer  is:  by  relaxation  iterations 
on  the  coarse  grid,  and  with  the  help  of  corrections  obtained  from  still  coarser  grids.  Thus,  in  a 
straightforward  recursive  fashion,  we  can  extend  the  above  two-level  equations  to  any  number  of 
levels.  In  view  of  (4)  and  the  fact  that  the  residual  for  the  level  k  equations  iB  given  by 


r**  =gh* -Afe*ufc\  (6) 

the  multi-level  equations  for  L  levels  are  given  by 

A^u**  =  %h *  for  1  <  k  <  L,  (7) 

where 


g'1'-  =  r*'-;  and 

gh‘  =  Aht(  I  ufc‘+>)+  I  (g'”"M  —  A'u+,u'“+1)>  for  0  <  fc  <  L  —  1 .  (8) 


Note  that  the  original,  right-hand  side  thk  of  the  level  problem  occurs  only  on  the  finest 
level  L  The  right-hand  sides  of  the  coarser  levels  have  been  modified  in  order  that  the  finest 
level  accuracy  be  properly  maintained  throughout;  that  is  to  say,  in  order  for  the  solutions  uhk 
to  coincide:  uhl  =  I*,***,  uha  =  •••  =  Ihj=»h,  •  •  -IkL  uh‘  -  Analogously  to  the  two-grid  case, 

we  can  interpret  the  difference  of  the  original  and  the  corrected  right-hand  sides,  thk  —  gh‘,  as 
an  estimate  of  the  local  truncation  error  of  level  it  relative  to  the  finest  level. 


5.3.  Multi-Level  Surface  Reconstruction  Algorithms 

We  have  motivated  the  multi-level  approach  to  surface  reconstruction  and  described  in  a 
quantitative  manner  its  basic  components  —  the  intra-level  relaxation  processes,  and  the  inter-level 
interpolation  processes.  It  now  remains  to  show  how  to  bring  the  components  together  into  an 
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algorithm  for  solving  the  multi-level  equations  (7)  and  (8).  Several  schemes  have  been  proposed 
[Brandt,  1977a,  1977b;  Brandt  and  Dinar,  1979j.  We  will  describe  one  which  is  appropriate  in  terms 
of  our  surface  approximation  problem.9  Before  defining  the  full  multi-level  surface  reconstruction 
algorithm,  we  will  define  its  main  procedure,  the  multi-level  surface  reconstruction  cycle. 

The  multi-level  surface  reconstruction  cycle  starts  at  the  (currently)  finest  level  l,  making 
several  cycles  to  the  coarser  levels  It  —  l  —  1,  l  —  2, . . . ,  1,  until  a  hierarchy  of  surface  representations 
which  are  as  accurate  as  is  possible  in  the  Shl  space  is  obtained.  Let  tk  denote  a  tolerance  for 
solving  the  equations  on  level  k.  £  and  r\  are  switching  parameters  which  are  given  appropriate 
values  below. 

Algorithm  I  —  Multi-Level  Surface  Reconstruction  Cycle 
Step  1  —  initialise  the  finest  level  l. 

Set  the  right  hand  side  of  the  level  1  problem  A.h‘uh‘  =  g*1'  to  the  original  right  hand 
side:  g,l‘  «—  I*1.  Introduce  the  initial  approximation  vhl  «—  Vq1  .  Set  £/  «—  0,'°  and  k  «—  l. 

Step  2  —  start  a  new  operation  level  k. 

Set  e£ld  «—  oo. 

Step  3  —  perform  a  relaxation  iteration. 

Perform  a  relaxation  iteration  for  the  equation  Ahk\thk  —  g*k  and  concurrently  compute 
some  norm  of  the  residual  given  by  (6),  e*  «—  ||rh,‘  ||. 

Step  4  —  test  the  convergence  and  its  rate. 

If  e*  <  tfc,  then  convergence  has  been  obtained  at  the  current  operation  level;  go 
to  Step  6.  If  k  =  1,  go  to  Step  3.  If  e*  <  r?e°ld  then  the  convergence  rate  is  still 
satisfactory,  set  e°'d  efc  and  go  to  Step  3;  otherwise  the  convergence  rate  is  slow  so 
go  to  Step  5. 

Step  5  —  transfer  to  coarser  level. 

Introduce  as  the  first  coarse-level  approximation  the  function  vhk~'  defined  by  the 
nodal  displacements 

, _  J  yAfc 

hk**hk—l 

Set  the  right-hand  side  of  the  coarser  level  problem  A.Kk-'uhk-'  =  gfck_1  to 
g'**-  —  A'^-'v'1*-*  +  I  (g'**  -  Ahkvhk) 


’Brandt  refers  to  it  as  the  accommodative,  full  multi-grid,  full  approximation  scheme  algorithm  [Brandt  and 
Dinar,  1979}. 

'“This  value  for  e<  is  temporary.  A  realistic  value  is  introduced  in  Step  5. 
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(in  view  of  equation  (8)).  Set  the  tolerance  c*~i  «-  Concurrently  with  the 
computation  of  gh<,~‘,  compute  the  norm  of  the  local  truncation  error  I**-1  —  g'u“* 
using  the  same  norm  as  in  Step  3.  If  k  —  I  set  ti  *-  —  gh“  'll  "  Finally,  set 

k  ♦-  k  —  1  and  go  to  Step  2. 

Step  6  —  use  converged  solution  uh“  to  correct  Oner  level. 

If  k  <  f,  make  the  correction  (in  view  of  equation  (5)) 

v^  +  i  *_  vA*+‘  -f-  I  (uh*  —  I  v/l*+1); 
set  k  «—  k  4  1  and  go  to  Step  2.  Otherwise,  if  k  -=  l  end. 


The  relaxation  operation  in  Step  3  can,  in  principle,  be  based  on  any  one  of  the  iterative 
methods  described  in  Appendix  B,  but  is  usually  a  simple  Jacobi  or  Gauss- Seidel  iteration.  For 
our  surface  reconstruction  problem,  in  view  of  equation  (4.9)  and  the  discussion  in  Appendix  B, 
the  Jacobi  relaxation  iteration  in  the  interior  of  0  is  given  by 


(<+!)_  _  1  f8/V'>  ■  v<*>  .-<*>  +v(>)  ^ 

*.j  -f  V'-1 ^  V,  +  l ^  v,o  — t  '  v>,3+i) 

-  +  vi’iliJ+l  4  v!'{liJ+1) 

^(v!-2,,  +  v!+2.,  +  v!',;-2  +  vl'.l  +  j)  + 


while  the  Gauss-Seidel  relaxation  iteration  is  given  by 


J'+l)  =  1 

&+0 


IfJ 

h 2 


(V'  +  ,)  4v<')  (v<,+,)4v{,)  ^ 

\v«  —  l,j  '  V«-H  •  v».j  —  1  ^  Y*,J  +  l  J 

_  ._<«)  .  v(>)  ) 
1  +  —  1  '  v*  —  1 , j -f- 1  '  v»  +  l  .J  +  iy 

+  ,  JO  ,„(>  +  ’)  ,  (.)  )  ,  1 

^2  \Vt  —  2,j  i  V«4-2  +  Vi.;4-2 )  »  » 


(9) 


(10) 


where  we  have  suppressed  the  superscripts  hi,  indicating  the  level,  and  have  instead  introduced  the 
bracketed  superscripts  which  indicate  the  iteration  number.  Analogous  formulas  for  the  boundary 
points  can  be  derived  from  the  computational  molecules  associated  with  the  boundary.  The  norm 
computed  in  Steps  3  and  5  can  be  the  discrete  Li  or  Lx  norm.  In  the  case  of  Gauss-Seidel 
relaxation,  it  is  quicker  to  compute  the  dynamic  norm,  as  the  iteration  progresses,  rather  than 
the  static  norm  (see,  e.g.  [Brandt,  1977b,  pg.  286]). 

Ar  important  feature  of  the  multi-level  algorithm  is  that  the  local  Fourier  analysis,  in 
addition  to  providing  a  prediction  of  the  convergence  rate,  enables  one  to  predict  near-optimal 
values  for  the  switching  parameters.  It  turns  out  that  the  actual  values  assigned  to  the  switching 
parameters  are  not  critical,  and  that  good  values  are  (  —  0.2  and  rj  =  p,  where  Ji  is  the  smoothing 

"The  constant  $  is  the  value  of  2  ~v  (see  (Brandt,  1979,  pg.  65]),  where  p  =  2  is  the  approximation 
order  of  the  secoud  partial  derivatives  of  the  energy  inner  product  that  is  achieved  by  our  quadratic 
element. 
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factor  of  the  relaxation  method  used  in  Step  3  of  the  algorithm  (see  (Brandt,  1977b,  pg  290]). 
The  order  of  the  interpolation  operators  is  determined  by  the  problem  itself;  i.e.  the  order  of 
derivatives  in  the  energy  inner  product.  For  the  coarse-to-fine  interpolation  =.»/,**,  in  Step 
6,  the  natural  second-degree  interpolation  of  the  element  polynomial  pk  may  be  used.  On  the 
other  hand,  simple  injections  perform  well  for  the  fine-to-coarse  transfers  I/, t  =»/,»_,  in  Step  5  and 
+  in  Step  6. 

Having  defined  the  multi-level  cycle  which  starts  at  the  finest  level  l  and  cycles  through  the 
coarser  levels,  we  will  employ  it  as  a  procedure  within  a  more  general,  full  multi-level  surface 
reconstruction  algorithm.  We  now  think  of  the  level  /  of  the  cycling  algorithm  as  the  currently 
finest  level ;  i.e.,  the  finest  level  for  which  an  approximate  solution  has  already  been  computed  by 
the  multi-level  cycle.  The  full  algorithm  works  in  the  opposite  direction,  the  currently  finest  level 
progressing  from’ the  coarsest  level  l  -=  1,  to  the  finest  level  l  =  L. 

Algorithm  2  —  Full,  Multi-Level  Surface  Reconstruction  Algorithm 

Step  1  —  solve  the  coarsest-grid  equation. 

Compute  by  relaxation  an  approximate  solution  u^’  to  the  coarsest-grid  equation 

\h'ah'  =  fh‘.  Set  l  «-  1. 

Step  2  —  set  a  new  finest  level  l. 

If  l  =  L  stop.  Otherwise  increment  the  currently  finest  level  l  l  -j-  1,  and  Bet  the  first 

approximation  on  the  new  level  to  be  the  function  Uq1  defined  by  nodal  displacements 

=  I*,.,-*,  uh,~l. 

Step  3  —  perform  a  multi-level  cycle. 

Invoke  Algorithm  1  and  when  it  ends,  go  to  Step  2. 

Note  that  the  solution  in  Step  1  will  be  performed  quickly  because  Sh'  has  relatively  few 
dimensions.  In  Step  3,  each  time  Algorithm  1  terminates  at  level  l,  we  have  obtained  a  hierarchy  of 
l  representations  whose  accuracy  is  the  best  possible  on  level  l.  The  currently  finest  approximation 
is  then  interpolated  to  the  next  finer  level  until  the  finest  level  L  is  reached.  Brandt  recommends 
a  somewhat  higher-order  interpolation  for  the  initial  interpolation  in  Step  2.  Third  order 

Lagrangian  interpolation  seems  adequate  for  our  surface  interpolation  problem,  as  we  will  see 
from  the  demonstrations  in  the  next  section. 

Algorithm  1  is  accommodative  in  that  it  makes  internal  checks,  based  on  the  computation 
of  norms,  to  determine  when  to  switch  levels.  For  many  types  of  problems,  accommodative 
algorithms  behave  in  a  fairly  fixed  manner,  performing  a  similar  number  of  iterations  on  each 
level  before  switching.  It  is  then  possible  to  avoid  computing  the  dynamic-residual  norm  in  Step 
3  of  Algorithm  1,  and  to  preassign  a  fixed  flow.  A  switch  to  the  coarser  level  Sh>—'  is  made  after 
ne  iterations  have  been  performed  at  level  5hV  Analogously,  a  switch  to  toe  finer  level  Shk+'  is 
made  after  n/  iterations  have  been  performed  on  level  Shk  since  the  last  return  from  the  finer 
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level.  nc  depends  or  the  smoothing  factor,  a  good  choice  being  nc  =  log.  1/  log  p.  Sometimes  nj 
varies  from  level  to  level.  For  a  more  extensive  discussion  see  [Brandt,  1979,  pg.  68  39].  Fixed 
algorithms  are  to  be  preferred  for  parallel  implementations  in  general,  and  from  a  biological  point 
of  view  in  particular. 

In  order  to  evaluate  the  performance  of  the  multi-level  surface  reconstruction  algorithm,  we 
define  a  unit  of  computation  called  a  work  unit  which  is  the  amount  of  computation  required  to 
perform  one  relaxation  iteration  on  the  finest  level  L.  It  is  roughly  equal  to  wNh'  ,  where  in  is 
the  number  of  operations  required  to  compute  the  residual  at  each  node12  and  Nh‘  is  the  number 
of  nodes  at  the  finest  level.  Since  there  are  one  quarter  as  many  nodes  on  level  k  —  1  as  there 
are  on  level  k,  only  1/4'  work  unit  is  required  to  perform  a  relaxation  iteration  on  level  L  —  t. 
The  proportionate  amount  of  computation  done  on  coarser  grids  thuB  diminishes  very  rapidly. 
Although,  for  accommodative  algorithms,  it  is  difficult  to  predict  the  total  number  of  operations 
consumed  by  the  inter-level  processes,  it  normally  turns  out  to  be  considerably  less  than  the  total 
inter-level  process  computation,  and  is  therefore  usually  ignored  (see  [Brandt,  1977a,  section  6.2]). 

A  final  issue  that  we  have  not  yet  considered  in  quantitative  terms  is  the  choice  of  appropriate 
values  far  the  (spring)  constant  0.  In  the  mathematical  and,  in  particular,  in  the  finite  element 
literatuie,  the  constraint  term  £5(t;)  of  equation  (2.2)  is  known  as  a  penalty  function  (see,  e.g., 
[Courant,  1943),  [Babuska,  1973),  (Strang  and  Fix,  1973),  [Zienkiewicz,  1974)).  The  incorporation 
of  penalty  functions  into  variational  principles  is  a  standard  way  of  approximately  satisfying 
essentia'  boundary  conditions  by  converting  them  into  appropriate  natural  boundary  conditions 
which  r.iay  be  handled  straightforwardly  by  the  finite  element  method.  Penalty  functions  are 
particularly  useful  when  the  essential  boundary  conditions  in  question  are  complicated,  or  when 
only  their  approximate  satisfaction  is  desired,  as  in  the  case  of  visual  surface  approximation.  An 
optimal  value  for  0  can  be  derived  through  the  following  considerations.  Let  w  be  the  solution 
to  our  surface  approximation  problem,  which  interpolates  the  known  depth  points.  As  usual, 
u  6  ^2(12)  denotes  the  exact  solution  to  the  variational  principle  (2.2),  including  the  penalty  terra 
(f5,  and  uh  €  Sh  denotes  the  finite  element  approximation  to  u.  Then,  there  will  be  a  balance 
between  the  error  w  —  ti  which  measures  how  closely  the  surface  fits  the  constraints  and  the  error 
n  —  uh,  due  to  minimizing  over  a  finite  element  subspace  [Strang  and  Fix,  1973,  pp.  132-133). 
Analyzing  this  balance,  Babuska  [1973)  determined  that  the  optimal  value  for  0  is  dependent  on  h 
and  is  given  by  0k  —■  ih~k,  where  7  is  a  constant  and  k  is  the  degree  of  the  complete  polynomial 
contain* d  in  Sh .  Therefoie,  for  our  quadratic  finite  elements,  k  —  2,  and  the  best  valu*  for  0  at 
level  j  of  the  multi-level  algorithm  is  0hl  —  "?/h2. 

5.4.  Examples  of  Multi-Level  Surface  Reconstruction 

Figure  7,  is  a  schematic  diagram  of  the  structure  of  the  multi-level  surface  reconstruction 

I2to  is  determined  by  the  specific  relaxation  scheme  used,  but  due  to  the  size  of  the  support  of  the 
central  computational  molecule,  it  is  approximately  equal  to  13. 
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Figure  7.  The  structure  of  the  multi-level  surface  reconstruction  algorithm. 


coarse 


medium 


fine 


TERZOPOULOS 


MULTI  LEVEL  SURFACE  RECONSTRUCTION 


algorithm,  showing  three  levels  of  resolution.  The  diagram  depicts  the  relaxation  processes  operating 
at  each  level,  as  well  as  the  fine-to-coarse  and  coarse-to-fine  processes  which  transfer  information 
between  levels.  The  algorithm  transforms  a  hierarchy  of  sparse  surface  depth  representations,  such 
as  might  be  provided  through  the  independent  stereo  bandpass  channels,  into  a  hierarchy  of  full 
surface  representations  which  constitute  the  full  2^-D  sketch.  The  constraints  for  the  surfaces 
shown  in  the  figure  are  random  samples  from  a  surface  which  varies  sinusoidally  in  depth.  It  is 
evident  that  the  multiple  full  representations  output  by  the  algorithm,  describe  the  sinusoidal 
surface  over  a  range  of  scales  and  that  the  accuracy  of  the  finest  representation  is  maintained  in 
the  coarser  ones. 

In  this  section,  a  number  of  examples  of  multi-level  surface  reconstruction  are  presented. 
We  will  consider  the  reconstruction  of  surfaces  from  artificially-generated  constraints,  as  well 
as  constraints  generated  from  natural  images  by  an  implementation  of  the  Marr-Poggio  stereo 
theory.  The  performance  of  the  multi-level  algorithm  is  compared  to  that  of  single-level  iterative 
algorithms.  In  the  examples  presented,  the  algorithm  was  started  from  identically  zero  initial 
approximations  on  all  the  levels. 

In  the  first  sequence  of  figures,  we  present  synthetic  examples  of  surface  reconstructions  with 
the  purpose  of  illustrating  the  performance  of  the  algorithm  in  reconstructing  quadric  surfaces 
having  :',ero,  positive,  and  negative  Gaussian  curvatures.  Constraints  on  each  level  were  synthesised 
by  sampling  depth  along  arcs  on  the  surface.  The  examples,  involved  four  levels  whose  grids  had 
dimensions  N*'  =  N$'  =  17,  N =  N =  33,  N =  65,  and  N**  =  Nhy'  =  129, 
with  corresponding  grid  spacings  hi  =  0.8,  ha  =  0.4,  /13  =  0.2,  and  h4  =  0  1.  The  relaxation 
method  employed  was  the  Gauss-Seidel  method  of  equation  (5.10),  and  a  value  of  2.0  vas  chosen 
for  7,  giving  phj  =  2.0/hJ. 

Figure  8  shows  depth  constraints  whose  values  are  consistent  with  a  cylindrical  surface  viewed 
at  four  resolutions.  The  constraints  lie  along  arcs  of  greatest  curvature.  Figure  9  illustrates  the 
hierarchy  of  full  surface  descriptions  reconstructed  by  the  four-level  algorithm.  Since  the  constraints 
on  all  the  levels  sample  the  same  ideal  cylindrical  surface,  the  full  surface  representations  coincide 
to  a  high  degree  of  accuracy.  Convergence  was  obtained  after  12.0  work  units.  For  comparison 
purposes,  the  finest  level  problem  was  isolated  from  the  coarser  levels  and  the  same  Gauss-Seidel 
relaxation  algorithm  was  applied  to  it.  Figure  10  shows  the  (single-level)  approximation  obtained 
after  800  work  units  (i.e.,  iterations).  It  is  clear  that  we  are  still  very  far  from  convergence. 
Although  the  approximation  is  generally  Bmooth,  it  has  large  low-frequcncy  error  components  and 
the  approximate  surface  lies  far  below  its  final  value  between  constraints  which  are  separated  by 
fairly  large  distances.  As  predicted  by  the  local  Fourier  analysis,  it  is  precisely  such  low  frequency 
error  components  that  local  iterative  algorithms  have  difficulty  liquidating.  In  fact,  the  following 
characteristic  phenomenon  was  observed.  During  the  initial  iterations,  the  corrections  made  to 
the  approximation  decreased  rapidly,  so  that  by  the  SOO1*1  iteration  they  are  minute,  even  though 
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Figure  10.  Single  level  approximation  after  800  work  units. 


the  error  norm  is  still  very  large.  Since  there  are  17, 361  nodes  in  the  grid,  it  may  take  on  the 
order  of  (17,361)2  work  units  to  obtain  the  solution  without  the  help  of  the  coarser  levels.  Thus, 
the  multi-level  algorithm  is  vastly  superior  when  the  constraints  are  far  apart. 

Figures  11  and  12  show  a  synthetic  example  of  the  reconstruction  of  a  hemispherical  surface 
from  constraints  which  form  latitudinal  circles.  The  hierarchy  of  full  surface  representations  was 
obtained  after  4.25  work  units.  Figures  13  and  14  illustrate  an  analogous  example  involving 
a  hyperbolic  paraboloid  (saddle  surface),  where  the  constraints  form  parallel  parabolic  arcs. 
Convergence  was  achieved  after  only  2.5  work  units  (only  a  single  iteration  was  performed  on 
the  finest  level).  Single  level  algorithms  applied  to  the  above  surfaces  exhibited  poor  convergence 
properties  similar  to  the  case  of  the  cylinder. 

The  above  examples  simulate  a  visual  situation  where  the  surface  in  the  scene  has  reflectance 
changes  in  the  form  of  widely- spaced  rulings  but  is  otherwise  free  of  intensity  changes.  This  is  an 
unlikely  situation  in  view  of  the  fact  that  the  visual  world  is  full  of  textures,  which  often  arise  from 
surface  material  and  pigment  changes.  Such  textures  generally  result  in  relatively  densely-spaced 
zero- crossings  forming,  to  a  certain  extent,  random  patterns.  In  turn,  these  xero- crossings  give  rise 
to  constraints  having  similar  properties.  Figure  15  illustrates  a  simulation  of  this  situation  using 
a  surface  varying  sinusoidally  in  depth.  A  three-level  surface  reconstruction  algorithm  was  used. 
The  constraints  input  on  each  level  were  30%-density  randomly-located  samples  of  the  surface 
depth.  In  addition,  to  simulate  the  effects  of  possible  inaccuracies  in  the  constraint  values,  each 
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Figure  15.  Surface  reconstruction  from  randomly-placed  depth  constraints. 


TERZOPOULOS 


MULTI -LEVEL  SURFACE  KKCONSTrX’CTION 


Figure  16.  Single  level  approximation  after  19  work  unit*. 


sample  was  corrupted  by  zero-mean,  uniformly-distributed,  additive  noise  whose  magnitude  was 
one  tenth  the  sample  value.  The  algorithm  generated  full  surface  representation  hierarchy  in  18.75 
work  units.  Evidently  our  spring  model  for  the  influence  of  the  constraints,  with  Pitj  =  2 /h2Jt 
is  adequate  for  this  case  in  that  the  additive  noise  has  not  adversely  affected  the  quality  of  the 
reconstructed  surfaces.  The  results  after  19  work  units  (iterations)  with  single-level  Gauss-Seidel 
relaxation  on  the  isolated  finest  grid  are  shown  in  Figure  16  for  comparison.  It  is  evident  that  the 
approximation  is  still  far  from  the  true  solution.  In  fact,  a  total  of  71  work  units  was  required 
to  reduce  the  error  norm  to  the  magnitude  obtained  after  18.75  work  units  by  the  three-level 
algorithm.  The  saving  in  computation  is  less  in  this  example  than  in  the  ones  above  because, 
first,  only  three  levels  were  used  and,  second,  the  density  of  the  constraints  is  greater.  In  general, 
the  greater  the  density  of  the  known  depth  values,  the  tighter  the  surface  is  constrained,  and 
the  convergence  is  expected  to  be  faster.  Another  way  to  think  of  this  is  that  as  the  average 
distance  between  constraints  decreases,  the  efficiency  of  relaxation  in  liquidating  the  low-frequency 
Fourier  components  in  the  error  increases  and,  therefore,  the  relative  advantage  of  the  multi-level 
algorithm  is  correspondingly  diminished. 

The  next  examples  illustrate  the  performance  of  the  multi-level  algorithm  using  disparity 
constraints  generated  by  Grimson's  implementation  of  the  Marr-Poggio  stereo  algorithm  [Grimson, 
1981a,  1981b],  which  includes  some  of  the  modifications  suggested  by  Mayhew  and  Frisby 
[1980,  1981]  for  exploiting  disparity  gradient  constraints  along  zero-crossing  contours.  The  stereo 
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I  algorithm  was  run  on  three  stereo  pairs  of  images,  shown  in  Figure  17,  which  were  digitized 

to  320  X  320  pixels  in  256  grey  levels.  The  pairs  from  top  to  bottom  are  A)  a  synthesized 
random  dot  stereogram  of  a  "wedding  cake”  of  stacked  planes,  B)  natural  images  of  a  portion 
of  a  coffee  jar  sprayed  with  “random  dot”  paint,  and  C)  an  aerial  view  of  a  set  of  buildings 
on  the  University  of  British  Columbia  campus.  A  three-channel  version  of  the  stereo  algorithm 
was  used.  The  resulting  sparse  disparity  representations  were  reduced  spatially  by  a  factor  of 
two,  and  input  to  a  three-level  surface  reconstruction  algorithm  whose  grids  had  dimensions 
7V£‘  =  =  41,  N*1  =  Ny1  —  81,  AfJ*  =  TVj3  =  161,  with  corresponding  grH  spacings 

hi  =  0.4,  h 2  =  0.2,  and  h.s  =  0.1.  The  surface  reconstructions  in  the  examples  is  based  on  the 
raw  disparities  whose  relation  to  depth  is  through  a  nonlinear  transformation.  Consequently  the 
-j  shapes  of  the  reconstructed  surfaces  are  distorted  to  a  certain  extent. 

,  The  sparse  disparity  constraints  provided  by  the  stereo  algorithm  and  the  hierarchy  of  full 

surface  representations  generated  by  the  three-level  reconstruction  algorithm  for  the  “wedding 
cake”  stereogram  are  shown  in  Figure  18.  The  value  7  =  0.5  was  used  in  the  algorithm,  and  the 
representations  shown  were  generated  after  16.75  work  units.  The  three-dimensional  structure  of 
the  planar  surfaces  is  clearly  evident  at  the  three  resolutions.  The  results  can  be  compared  to 
j  Figure  19  which  shows  the  approximation  obtained  by  a  single-grid  algorithm  on  the  finest  grid, 

after  more  than  900  work  units. 

Figure  20  illustrai  s  the  sparse  constraints  and  the  full  surface  representations  obtained 
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Figure  18.  Disparity  constraints  ami  full  surface  representations  of  “wedding  cake' 
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Figure  21.  Disparity  constraints  and  surfaces  reconstructed  Irom  the  aerial  view. 
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(with  7  =  0.1)  from  the  images  of  the  coflce  jar.13  The  reconstruction  required  1G.0  work  units. 
Finally,  Figure  21  shows  the  sparse  constraints  and  the  reconstruction  obtaiued  from  the  aerial 
view  (with  7  =  8.0),  after  21.875  work  units.  The  representations  are  displayed  as  grey  level 
images,  in  which  darkness  is  proportional  to  disparity. 

It  should  be  noted  from  the  above  examples  that  the  multi-level  surface  reconstruction 
algorithm,  in  its  present  form,  attempts  to  reconstruct  a  single  surface  over  the  entire  grid.  As 
a  consequence,  serious  problems  arise  near  sharp  changes  in  depth  such  as  those  due  to  partial 
occlusions  of  surfaces  in  the  scene.  The  reconstructed  surface  gives  the  undesirable  impression 
of  a  tablecloth  thrown  over  a  3-D  model  of  the  scene.  The  source  of  this  difficulty  is  discussed 
further  in  the  concluding  chapter,  where  possible  ways  of  overcoming  it  are  suggested. 

6.  GENERALIZED  INTERPOLATION  PROBLEMS  IN  VISION 

The  key  to  the  solution  of  many  problems  in  early  vision  is  the  imposition  of  constraints  based 
upon  assumptions  about  the  visual  world  which  are  almost  always  true.  A  common  assumption  is 
that  matter  is  cohesive;  i.e,  that  surfaces  are  continuous  over  most  of  the  Beene.  This  assumption 
is  usually  introduced  in  the  form  of  smoothness  constraints,  such  as  those  characterizing  the 
most  consistent  surface  in  visual  surface  reconstruction.  From  our  study  of  this  problem,  we 
have  seen  that  it  is  appropriate  to  formulate  smoothness  constraints  within  variational  principles. 
In  this  chapter,  we  study  a  general  class  of  variational  principles,  and  we  propose  that  the 
functionals  characterizing  these  variational  principles  are  appropriate  semi-norm b  for  formulating 
smoothness  constraints  because  they  possess  several  invariance  properties  which  become  important 
in  applications  to  early  vision.  In  order  to  simplify  the  discussion,  the  analysis  will  be  in  terms  of 
interpolatory  constraints  and  domains  of  infinite  extent.  Our  visual  surface  reconstruction  problem 
will  be  shown  to  be  a  special  case  of  this  generalized,  optimal  interpolation  problem  which  is  a 
natural  generalization  of  the  familiar  curve-fitting  problem  involving  splines. 

The  classical  spline  problem  involves  the  minimization  of  the  quadratic  functional 

|v 

under  the  interpolatory  constraints  v(x,)  =  c,,  1  <  «  <  Nc  with  Ne  >  m  >  1,  where  the  s, 
are  given  distinct  points  in  [a,  6]  and  the  c,  are  given  real  scalars.  The  natural  setting  for  this 
problem  is  a  vector  space  V  formed  by  the  class  of  functions  whose  (distributional)  derivatives 
up  to  order  m  are  in  Lj(a,  6);  that  is,  the  class  of  functions  which  are  elements  of  the  Sobolev 
space  of  order  m  over  [a,  6],  l/m([o,  6]),  defined  in  Section  A.l.  |  |m  is  a  semi-norm  which  is 
derived  from  a  semi-inner  product  and,  equipped  with  it,  V  becomes  a  semi-Hilbert  space.  The 

l3In  this  example,  the  constraints  for  the  coarser  channels  were  generated  by  averaging  down  the 
finest-channel  disparities. 


-/ 
•  a 


dTv^x) 


dx " 


dx 


(1) 
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conditiooi  imposed  on  the  constraints  ensure  the  existence  of  a  unique  solution  it  Q  S,  where  S 
is  the  convex  subset  of  V  whose  elements  interpolate  the  constraints.  The  characterization  of  u 
as  an  odd-degree  polynomial  spline  and  various  intrinsic  properties  such  as  the  minimum  norm 
and  best  approximation  properties  follow  from  the  orthogonal  projection  theorem  [Ahlberg  et  al., 
1967]  (see  also  the  proof  of  Theorem  A.l). 

Duchon  [1976,  1977]  and  Meinguet  [1979a,  1979b]  describe  an  n-dimensional  generalization 
of  the  optimal,  univariate  spline  interpolation  problem.  The  generalized  optimal  interpolation 
problem  involves  the  minimization  of  the  functional  j-|?n,  where 


-(.  t  ,/, 


a’“u(x) 

dxti  . . .  dx,„ 


'■  \ * 
dxj  , 


(2) 


and  m  and  n  arfe  given  positive  integers.  The  generalized  interpolation  problem  is  naturally  set 
in  a  space  V  of  functions  which  are  elements  of  the  Sobolev  space  jf’"(!Rn).14  |  |m  is  a  semi-norm 
whose  mi //-space15  is  the  M  —  (n+™"~1)  dimensional  space  of  all  polynomials  over  9?"  of  degree 
less  than  or  equal  to  m  —  1:  M  =  IIm— '(SR”)  [Schwartz,  1966,  pg.  60].  Equipped  with  the 
semi-inner  product  corresponding  to  |-|m,  V  becomes  a  soini-Hilbert  space.16 

Let  the  finite  set  of  distinct  constraints  C  ~  {(x,,c,)  |  1  <  t  <  Nc,  x,  £  SR",  c,  £  3?}  contain 
a  subset  of  M  members  such  that  there  exists  a  unique  element  p£  IT"-1  (SR")  of  the  null  space 
of  |  |m  which  interpolates  the  M  constraints  in  the  subset;  that  is,  such  that  there  exists  a  unique 
polynomial  of  degiee  m  —  1  which  satisfies  the  conditions  p(ij)  =  c;  for  each  j  which  indexes  a 
constraint  of  the  above  subset.  We  call  such  a  subset  an  M -un’rolvent  subset. 


We  can  pose  the  generalized  optimal  interpolation  problem  in  the  following  way.  Given  a 
set  of  constraints  C  which  contain  an  A/- unisolvent  subset,  find  that  element  u  £  S  such  that 

Mm  =  inf  Ml. 

v£S 

where  once  again  S  is  the  set  of  functions  in  V  which  interpolate  the  constraints.  The  problem  is 
well-posed  because  we  are  minimizing  a  semi-norm  within  a  convex  subset  S  of  a  semi-Hilbert  space 
and,  furthermore,  the  existence  of  an  .V-unisolvent  subset  of  constraints  reduces  the  null-space  of 
the  functional  to  at  most  a  single  nonzero  element  of  5.  A  solution  u  is  then  guaranteed  to  exist 
and  be  unique  by  the  orthogonal  projection  theorem  (see  the  proof  of  theorem  A.l).17 

l4More  precisely,  the  Bpace  V  is  the  Beppo  Levi  space  [Duchon,  1977;  Meinguet,  1979a,  1979b]  of  order 
m  over  !?n  defined  by  BLm(Sn)  =  {v  |  £  Lj  for|a|  =  m},  where  a  =  (ai,...,a„)  is  a  “vector"  of 

positive  ntegers  and  |o|  =  oi  -f-  ■  •  ■  -f-  a,;  that  is,  it  is  the  vector  space  of  functions  for  which  all  the 
partial  derivatives  of  (total)  order  m  are  square  integrable  in  !Rn.  The  Beppo  Levi  spaces  are  related  to 
the  Sobolev  Spaces. 

r’The  cull  space  of  a  semi-norm  is  the  space  of  functions  which  the  semi-norm  maps  to  sero. 

''According  to  the  Sobolev  inequality  (see  section  A.l),  when  m  >  n/2,  V  is  a  semi-Hilbert  function 
space  of  continuous  functions  (Meinguet,  1979a,  1979b). 

17Moreover,  as  a  consequence  of  the  Sobolev  inequality  given  in  Section  A.l,  is  will  be  continuous  if 
m  >  n/2. 
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Why  is  the  class  of  semi-norms  defined  in  (2)  important  in  the  context  of  vision?  As  was 
argued  recently  by  Brady  and  Horn  [1981],  many  processes  in  early  vision  are  approximately 
isotropic  and,  therefore,  it  seems  that  operators  which  model  these  processes  ought  to  be 
rotationally  symmetric.  An  example  of  such  an  operator  is  the  V 2 G  edge  operator  proposed  for 
computing  the  primal  sketch  [Marr  and  Hildreth,  1980;  Hildreth,  1980].  The  class  of  semi-norms 
defined  in  (2)  is  of  interest,  since  all  its  members  |uj,„  are  invariant  under  rotation  and  translation 
transformations  and,  moreover,  if  a  dilation  or  contraction  in  \s  is  applied  to  v,  they  are 
multiplied  by  some  power  of  |\|  [Duchon,  1977,  pg.  86].  Therefore,  corresponding  interpolation 
methods  will  commute  with  any  similarity  transformations  applied  to  the  constraints.  Clearly, 
these  properties  arc  essential  for  interpolation  processes  which  contribute  to  th2  generation  of  the 
2^-D  sketch  —  the  surfaces  generated  by  surface  reconstruction  algorithms  should  not  change 
shape  as  the  objects  in  the  scene  undergo  translations  or  rotations  parallel  to  Ihe  image  plane,  or 
undergo  displacements  directly  towards  or  away  from  the  viewer. 


For  certain  instances  of  m  >  1  and  n  >  1,  the  general  interpolation  problem  has  familiar 
physical  interpretations  which  are  most  often  encountered  in  a  differential  form  through  the 
associated  Euler-Lagrange  equations.  Consider  first  the  one-dimensional  case,  n  =  1.  The 
generalized  interpolation  problem  then  reduces  to  the  common  univariate  spline  problem  of 
equation  (1).  The  particular  value  chosen  for  m  determines  the  order  of  continuity  of  the  optimal 
curve  —  as  m  increases,  the  smoothness  of  the  solutions  increases.  In  particular,  for  the  case 
m  ==  1,  |u|J  =  J,x  v\  dx  measures  the  ene;r„>  n  a  string  of  infinite  extent,  and  leads  to  interpolants 
having  C°  continuity.  The  associated  Euler-i.agrang  equation  is  uXI  =  0.  C'  continuity  may 
be  imposed  on  the  interpolant  by  choosing  m  2.  In  this  case,  ]u ||  =  v^7  dx  measures  the 

strain  energy  of  bending  in  a  thin  beam  of  infinite  extent,  and  the  Euler-Lagrange  equation  is 
uIZZx  =  0.  This  class  of  univariate  semi-norms  seems  to  be  appropriate  for  imposing  continuity 
constraints  in  the  computation  of  optical  flow  along  zero-crossing  contours  in  the  primal  sketch 
[E.C.  Hildreth,  personal  communication]. 


Next  consider  the  generalized  interpolation  problem  in  two  dimensions.  For  n  =  2,  the 
semi-norms  become 


m,  once  again,  determining  the  degree  of  smoothness  of  the  solution.  For  m  —  1,  | vjf  = 
//■•(  vj  +  v2)  dx  dy  measures  the  potential  energy  related  to  the  statics  of  a  membrane  (rubber 
sheet),  and  the  associated  Euler-Lagrange  can  be  shown  to  be  Laplace’s  equation,  Au  =  0 
[Courant  and  Hilbert,  1953,  pg.  247].  A  semi-norm  of  this  order  implicitly  imposes  the  smoothness 
constraints  in  algorithms  proposed  for  computing  lightness  [Horn,  1974],  shape  from  shading 
[Ikeuchi  and  Horn,  1981],  optical  flow  [Horn  and  Schunck,  1981],  photometric  stereo  [lkeuchi, 
1981],  etc.  With  m  =  2,  the  smoothness  of  the  interpolating  surface  is  increased  to  C1,  the 
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functional  taking  the  form  |t»|^  =  /  +  2v“v  try)dzdy.  This  will  be  recognized  as 

being  our  familiar  functional  representing  the  strain  energy  of  the  thin  plate  £1  with  Poisson 
constant  a  =  0  (refer  to  equation  (2.1)),  whose  Euler-Lagrange  equation  was  shown  to  be  the 
biharmonic  equation,  A2u  =  0.18  As  we  have  demonstrated,  this  order  of  smoothness  seems  to 
be  most  appropriate  in  visual  surface  reconstruction  from,  e.g.,  stereo  information  (refer  also  to 
the  discussion  on  the  “quadratic  variation’’  in  {Grimson,  1981a]). 

It  becomes  clear  that  we  are  dealing  with  a  class  of  quadratic  variational  problems,  the  order 
of  whose  Euler  equations  is  determined  by  the  degree  of  smoothness  demanded  of  the  solutions.  For 
m  =  1  we  obtain  Laplace’s  equation  in  n  dimensions,  for  m  =  2  the  biharmonic  equation,  and  so 
on.  In  general,  the  Euler  equation  is  an  n-dimensional,  linear,  elliptic  partial  differential  equation 
of  order  2m.  Moreover,  the  general  interpolation  problems  have  straightforward  formulations 
as  analogous  approximation  problems.  For  example,  we  can  define  appropriate  constraint  terms 
analogous  to  the  term  f5  (equation  (2.2))  for  our  surface  approximation  problem.  Hence,  there 
exists  a  general  framework  in  which  to  solve  functional  approximation  problems,  of  the  type 
arising  naturally  when  imposing  smoothness  constraints  in  early  vision.  Meaningful,  problems  can 
be  formulated  in  any  number  of  dimensions,  and  the  degree  of  smoothness  that  the  solutions 
should  possess  can  be  specified  o  priori.  In  this  sense  then,  the  Sobolev  spaces  can  be  viewed 
as  ingenious  formalizations  of  the  notion  of  the  “degree  of  smoothness”  of  admissible  functions 
and  therefore  are  ideal  domains  in  which  to  pose  and  solve  these  problems.  By  specifying  the 
(order  of)  the  Sobolev  space  to  which  the  solution  should  belong,  we  designate  its  position 
in  the  wide  spectrum  from  very  smooth  functions  to  singular  distributions.  Satisfaction  of  the 
requirements,  that  the  admissible  space  be  a  semi  Hilbert  space  and  that  the  constraints  include 
an  .A/-unisolvent  subset  will  guarantee  uniqueness.  Needless  to  say,  the  theory  of  the  finite  element 
method  is  applicable  to  either  the  interpolation  or  approximation  formulations,  and  it  will  dictate 
appropriate  finite  element  discretization  schemes  for  the  associated  variational  principles. 

When  solving  these  variational  principles  using  local,  iterative  algorithms  such  as  the 
ones  described  in  this  paper,  smoothness  constraints  are  imposed  globally  over  retinocentric 
representations  by  a  process  of  constraint  propagation.  Inspired  by  the  work  of  Waltz  [1975],  a 
class  of  algorithms  called  relaxation  labeling  algorithms  were  introduced  as  cooperative,  constraint 
propagation  techniques  in  vision  and  image  processing  by  Rosenfeld,  Hummel,  and  Zucker  [1976]. 
Although  they  have  seen  extensive  use  [Davis  and  Rosenfeld,  1981],  their  generality  has  made 
them  difficult  to  understand  and,  unlike  the  techniques  and  algorithms  which  are  the  subject  of 
this  paper,  the  foundations  of  most  relaxation  labeling  schemes  are  unfortunately  poorly- developed 
mathematically. 


'*Duchon  [1977]  understandably  refers  to  the  solutions  as  thin  plate  splines  which  also  reflects  the 
fact  that  they  are  natural  two-dimensional  generalizations  of  commonly-used  univariate  splines.  In  the 
engineering  literature  .they  are  called  surface  splines  [Harder  and  Desmarais,  1972]. 
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Recently,  some  theoretical  understanding  has  been  achieved  by  viewing  relaxation  algorithms 
as  techniques  for  solving  constrained  optimization  problems  (see,  e.g.,  [Ullman,  1979b],  [Faugcras 
and  Bcrthod,  1981],  [Hummel  and  Zucker,  1980]).  From  this  new  point  of  view,  the  relationships 
between  relaxation  labeling  techniques  and  iterative  solution  of  finite  element  equations  arising 
from  variational  formulations  become  clearer  —  relaxation  labeling  schemes  can  be  viewed  as 
iterative  algorithms  for  solving  optimal  approximation  problems  over  closed  convex  subsets  (of 
possible  labelings)  [Hummel  and  Zucker,  1980],  Necessary  conditions  for  solutions  (fixed  points) 
are  then  expressed  as  sets  of  variational  inequalities  [Ciarlet,  1978;  Kinderlehrer  and  Stampacchia, 
1980]  and  appropriate  updating  rules  are  natural  generalizations  of  the  classical  local  iterative 
methods  for  solving  large  systems  of  linear  [Young,  1971]  or  nonlinear  [Ortega  and  Rhcinboldt, 
1970]  equations.  Moreover,  if  the  compatibility  functions  among  neighboring  nodes  are  symmetric, 
then  there  exist  associated  variational  principles  defining  equivalent  formulations  as  minimization 
problems.  Fortunately,  it  is  possible  to  apply  the  finite  element  method  to  nonlinear  problems 
stemming  from  variational  inequalities  [Ciarlet,  1978].  In  a  certain  sense  then,  finite  elements  can  be 
viewed  as  systematically-derived,  physically-based  compatibility  relationships  among  neighboring 
nodes.  In  view  of  the  relationships  between  the  two  techniques,  it  is  hoped  that  aspects  of  our 
multi-level  approach  to  solving  the  discrete  finite  element  equations  for  the  surface  reconstruction 
problem  may  contribute  to  the  theory  of  hierarchical  relaxation  labeling  [Davis  and  Rosenfeld, 
1 981;  Zucker,  1978;  Zucker  and  Mohammed,  1979]. 

7.  SUMMARY  AND  EXTENSIONS 

Information  about  the  shapes  of  visual  surfaces  that  is  inferred  from  the  retinal  images  in 
the  early  computational  stages  in  vision  is  sparse.  Nevertheless,  our  perception  is  that  of  full 
(piecewise)  continuous  surfaces.  In  this  paper  we  have  proposed  a  hierarchical  approach  to  the 
reconstruction  of  full  surface  representations,  consistent  with  our  perception  of  the  visual  world. 
The  foundations  of  our  paradigm  are  embedded  in  a  tight  mathematical  formalism  which  at 
the  same  time  seems  sufficiently  general  to  encompass  many  aspects  of  the  complex  information 
processing  task  which  is  the  generation  of  the  full  2^-D  sketch. 

Visual  surface  reconstruction  was  formulated  as  an  optimal  approximation  problem  having 
an  intuitively  simple  physical  interpretation  —  a  thin  flexible  plate  which  is  allowed  to  achieve  sn 
energy-minimizing  state  of  stable  equilibrium  under  the  influence  of  externally-imposed  constraints. 
This  physical  model  led  directly  to  an  analysis  in  terms  of  the  calculus  of  variations,  and  a 
proof  that  the  problem  is  well-posed  in  practice.  The  model  also  suggested  a  class  of  techniques 
for  optimally  approximating  the  continuous  solution  by  an  equivalent  discrete  problem  which  is 
amenable  to  computational  solution.  We  chose  to  apply  the  finite  element  method  for  reasons 
which  include  its  generality,  the  availability  of  a  tight  theory  governing  its  use,  the  simple 
discrete  problem  to  which  it  gives  rise,  and  its  promise  in  vision  as  a  systematic  methodology 
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for  constructing  local  representations  of  surfaces. 19  At  each  step,  the  underlying  mathematical 
theory  assured  us  that,  ultimately,  our  problem  would  have  a  unique  solution  that,  in  principle, 
could  be  computed  by  biologically- based  mechanisms.  Our  search  for  efficient  algorithms  and 
our  insights  into  the  multi-level  structure  of  the  early  processing  stages  in  vision  led  us  to  a 
multi-level  algorithm  which  solves  simultaneously  a  hierarchy  of  surface  approximation  problems 
spanning  a  range  of  resolutions.  The  local-support  processes  comprising  the  algorithm  include 
iterative  intra- level  relaxation  processes,  and  inter-level  processes  which  serve  to  communicate 
information  between  levels.  The  inter-level  process  include  injections  from  (ine  grids  to  coarse 
and  polynomial  interpolations  from  coarse  grids  to  fine.  Tests  on  stereo  data  verified  that  our 
multi-level  surface  reconstruction  algorithm  meets  theoretical  expectations  of  increased  speed  and, 
moreover,  generates  a  potentially  useful  hierarchy  of  consistent  surface  representations.  Finally, 
we  examined  our  basic  surface  approximation  problem  in  a  more  general  setting,  and  related  it 
to  a  broader  class  of  optimal  approximation  problems  based  on  semi-norms  that  commute  with 
similarity  transformations  applied  to  the  constraints,  a  property  which  is  important  in  the  context 
of  vision. 

Although  we  have  laid  down  the  foundations  of  our  approach  primarily  in  terms  of  stereopsis, 
the  methodology  is  by  no  means  limited  to  the  type  of  information  produced  by  this  particular 
module.  Indeed,  our  point  of  view  speaks  to  the  broader  issue  of  how  to  combine  the  information 
about  the  shapes  of  visible  surfaces  generated  by  various  vision  modules  into  a  self-consistent 
whole.  Several  possibilities  arise,  some  of  which  we  will  now  consider  briefly. 

The  simultaneous  assimilation  of  information  from  different  sources  can  be  realized  by 
defining  more  sophisticated  penalty  functions  to  replace  £5  in  Chapter  2.  For  example,  in  the  case 
of  depth  constraints  from,  say,  stereo  and  motion,  we  can  straightforwardly  introduce  additional 
terms  of  the  same  form  as  £5  for  each  process.  In  terms  of  our  plate  model,  we  introduce  two  sets 
of  imaginary  pins  with  attached  springs,  and  allow  the  possibility  of  a  constraint  generated  by 
stereo  to  coincide  with  one  provided  by  motion.  Imperfections  in  the  retinal  images  are  likely  to 
affect  the  two  processes  in  different  ways,  and  moreover  each  will  in  its  own  way  sporadically  fall 
prey  to  gross  misinterpretations  of  the  information  in  the  primal  sketches.  Whatever  the  situation, 
our  physical  model  assumes  an  energy-minimizing  state,  and  the  resulting  surface  is  an  optimal 
compromise  in  view  of  the  constraints  provided.  In  places  where  the  information  is  consistent,  the 
final  interpretation  is  reinforced.  In  places  where  there  is  a  conflict,  it  is  resolved  by  competition 
with  nearby  constraints  from  both  processes. 

The  influence  of  each  constraint  may  be  controlled,  possibly  dynamically  by  the  processes 
themselves,  by  assigning  different  values  to  each  spring  constant.  For  example,  different  confidence 
value »  may  be  given  to  individual  constraints  generated  by  the  stereo  matcher,  according  to 

’"On  (fail  latter  point  we  should  mention  again  that  the  finite  element  method  allows  us  to  handle 
domains  of  complex  shape,  natural  boundary  conditions,  and  to  set  up  nonuniforra  discretizations  of  the 
domain  —  e.g.,  to  vary  the  resolution  across  the  domain. 
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regional  statistics  of  the  rate  of  successfully  matched  zero-crossings,  which  may  haw;  to  be 
computed  anyway  [Marr  and  Poggio,  1979;  Grimson,  1981a,  sections  2.5  and  3.1).  The  extent  of 
the  constraint's  influence  on  the  surface  may  also  be  varied  by  extending  our  model  to  that  of 
an  inhomogeneous  plate,  whose  flexibility  varies  over  the  domain.  Numerous  possibilities  exist  for 
defining  weighting  functions  to  apply  to  the  strain  energy  density  of  the  plate.  All  such  proposals 
for  modifying  the  form  of  the  functional  must  be  first  be  shown  to  lead  to  well-posed  problems, 
by  extensions  of  the  analysis  carried  out  in  Chapters  3  and  4. 

Another  important  issue  is  how  to  incorporate  information  other  than  measurements  of  the 
distance  to  the  surface.  An  important  class  of  processes  generate  cues  about  visual  surfaces  in 
the  form  of  local  orientation  measurements.  Examples  in  this  class  are  the  analysis  of  occluding 
contours  [Marr,  1977],  as  well  as  “shape-from”  processes  such  as  shape  from  shading  [Horn,  1975], 
contours  and  texture  [Kender,  1980;  Stevens,  1981;  Witkin,  1981],  regular  patterns  [Kanade, 
1981],  etc.  The  finite  element  method  provides  a  general  way  of  handling  orientation  constraints 
through  the  use  of  elements  with  degrees  of  freedom  which  include  the  first  partial  derivatives  of 
the  surface.  An  appealing  example  is  Adini’s  rectangle  which  was  described  in  Section  4.1.  Surface 
representations  based  on  this  element  would  make  explicit  the  information  about  the  local  slope 
of  surfaces,  as  well  as  their  distance  from  the  viewer.  Discrete  problems  derived  by  applying  this 
element  would  correspond  to  a  coupled  system  of  two  discrete  Euler-Lagrange  equations  for  the 
plate,  a  fourth-order  equation  for  the  displacements  at  the  nodes,  and  a  second-order  equation 
for  the  slopes.  On  the  other  hand,  we  pay  a  price  for  this  added  capability  —  the  dimension  of 
the  finite  element  space  is  tripled,  making  the  resulting  discrete  system  even  larger.  Nevertheless, 
the  price  may  turn  out  to  be  worth  paying  in  order  to  obtain  useful  surface  representations 
and,  moreover,  it  may  not  be  too  high  when  massively  parallel  computational  mechanisms  are 
contemplated. 

A  different  way  of  handling  orientation  constraints  which  can  be  used  with  our  simple 
quadratic  elements  is,  once  again,  by  the  use  of  appropriate  penalty  functions.  In  terms  of  our 
model,  we  can  imagine  the  situation  for  a  single  constraint  and  a  surface  patch  as  illustrated 
in  Figure  22.  Here,  we  attach  a  spring  between  the  surface  normal,  an  imaginary  quill  rigidly 
fixed  to  the  plate’s  surface  at  a  particular  point,  and  the  orientation  constraint,  another  quill 
emanating  from  the  same  point,  but  having  a  fixed  orientation  in  space.  Given  this  arrangement, 
the  surface  is  “pulled”  locally  so  that  its  orientation  tends  to  align  with  that  of  the  constraint. 
The  appropriate  penalty  term  is  the  potential  energy  in  the  spring.  This  energy  can  be  expressed 
straightforwardly  in  terms  of  the  first  partial  derivatives  of  the  quadratic  surface  patch  within 
the  element,  and  ultimately  in  terms  of  the  node  displacements. 

A  more  immediately  important  issue,  one  that  was  raised  in  reference  to  the  examples  of 
surface  reconstruction  presented  earlier,  is  that  of  dealing  with  depth  die  continuities.  In  its  present 
form,  the  surface  approximation  algorithm  can  deal  in  a  meaningful  way  with  scenes  containing 
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Figure  22.  Physical  model  of  the  effect  of  an  orientation  constraint. 


orientation  constraint 


only  a  single  surface.  This  is  due  to  the  fact  that  it  does  not  incorporate  the  notion  of  an  occluding 
contour ;  that  is  to  say,  it  attempts  to  fit  a  single  surface  over  the  whole  sparse  2)-D  sketch, 
interpolating  indiscriminately  across  contours  which  correspond  to  places  where  surfaces  in  the 
scene  occlude  one  other  from  the  viewer.  Clearly,  this  action  is  inappropriate  since  the  surfaces 
on  either  side  of  the  occluding  boundary  ought  to  have  no  influence  on  one  another.  Moreover, 
the  variational  principle  for  surface  approximation  was  based  on  the  small  deflection  theory  of 
the  plate20  and,  consequently,  we  expect  our  surface  to  behave  strangely  in  the  vicinity  of  a  large 
change  in  depth,  resulting  in,  for  example,  a  Gibb’s  phenomenon  similar  to  that  observed  when 
approximating  discontinuous  functions  with  Fourier  series. 

How  can  these  depth  discontinuities  be  detected  and  how  do  we  prevent  interpolation 
across  them?  Grimson  [1981a,  section  9.4]  noted  the  importance  of  this  question  and  made  some 
speculations  about  possible  answers  to  the  first  of  its  two  parts.  The  feasibility  of  his  suggestions 
remains  an  open  question.  Here,  we  would  like  to  propose  another  approach  that  is  again  suggested 
by  our  physical  model.  In  places  of  sharp  changes  in  depth  (or  surface  orientation)  the  strain 
energy  in  the  plate  will  be  locally  high.  Measuring  this  energy  locally  is  a  simple  matter  —  we  use 
the  energy  inner  product  to  compute  the  strain  energy  norm  os( vh,  vk)$  over  the  element  domains. 
Points  of  high  strain  energy  are  likely  candidates  Tot  inferring  the  presence  of  discontinuities  in 

30Tbe  large  deflection  plate  bending  theory  is  considerably  more  complicated  and  leads  to  an 
Eulcr-Lagrange  equation  in  the  form  of  two  coupled,  nonlinear,  fourth-order  partial  different  al  equations 
known  as  von  Karmann’s  equations  (see,  e.g.,  [Landau  and  Lifshiti,  1970)  or  [Mansfield,  1964)) 
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depth.  We  can  also  exploit  our  expectations  about  the  world  for  added  constraint,  and  assert 
that,  since  most  of  the  retinal  image  is  made  up  of  coherent  surfaces,  occlusions  in  depth  are 
likely  to  form  contours  in  the  image  and  not  be  sparsely-occurring  points.  Hence,  we  look  for 
contours  along  the  surface  of  the  plate,  where  the  strain  energy  is  high.  Having  located  the 
occluding  contours,  the  answer  to  the  second  part  of  the  question  is  simple,  in  principle.  To 
prevent  interpolation  across  different  surfaces,  we  “break”  the  plate  along  occluding  contours. 
Mathematically  speaking,  this  is  done  by  removing  plate  elements  along  the  occluding  contour, 
thereby  introducing  free  boundaries. 

Our  multi-level  approach  to  surface  reconstruction  constitutes  a  computational  paradigm 
which  has  contributed  toward  a  more  complete  understanding  of  the  generation  of  the  2£-D  sketch. 
Many  details  such  as  the  combination  of  information  from  the  various  modules  in  early  vision  and 
the  isolation  of  depth  discontinuities  remain  to  be  worked  out  rigorously  within  the  paradigm.  In 
addition,  a  number  of  exciting  issues  are  raised.  For  example,  how  can  the  hierarchy  of  surface 
representations  generated  by  the  algorithm  be  used  to  advantage  during  later  computational 
stages  in  which  three-dimensional,  object-centered  representations  are  generated  and  objects  are 
recognized.  Similar  implications  directed  to  the  related  field  of  robotics  and  manipulation  atao 
suggest  themselves.  Research  addressing  some  of  these  issues  is  currently  in  progress. 
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A.  THE  FINITE  ELEMENT  METHOD 


When  it. is  impossible  to  derive  an  analytical  solution  to  a  continuous  variational  principle, 
it  is  usual  to  attempt  an  approximation  by  defining  a  discrete  problem  which  is  similar  to  the 
continuous  one  and  which  leads  to  a  discrete  solution.  To  this  end,  we  will  first  state  an  abstract 
variational  principle  which  will  lead  us  to  an  optimal  approximation  to  the  exact  solution.  The 
variational  principle  is  called  abstract  inasmuch  as  it  represents  a  formulation  which  is  common  to 
a  variety  of  physical  problems,  such  a;  the  physical  model  for  our  surface  reconstruction  problem. 
We  will  also  state  theorems  which  give  conditions  guaranteeing  the  existence  and  uniqueness  of 
the  approximate  solution  and,  in  addition,  we  will  discuss  the  optimal  properties  of  the  proposed 
approximation. 


The  abstract  variational  principle  and  the  associated  theorems  are  stated  in  a  form  which  is 
convenient  for  the  application  of  the  finite  element  method,  a  powerful  technique  for  obtaining,  by 
numerical  means,  discrete  solutions  to  variational  problems.1  The  following  sections  develop  the 
mathematical  machinery  which  we  will  require  to  successfully  apply  the  method.  Key  mathematical 
ideas  include  a  set  of  Hilbert  spaces  (th^  Sobolev  spaces)  and  their  norms,  a  bilinear  form  (the 
energy  inner  product)  which  is  naturally  associated  with  the  specific  problem,  and  certain  optimal 
properties  of  the  (Ritz)  approximation  over  finite  dimensional  subspaces.  These  ideas  lead  to 
a  clean  and  precise  theory  governing  the  application  of  the  finite  element  method,  even  for 
complicated  geometries.  Comparatively  tight  theories  are  unavailable  for  alternate  approximation 
techniques  which  naturally  arise  from  nonvariatioria)  problem  statements;  e.g.,  the  finite  difference 
method  which  can  be  applied  to  equivalent  formulations  in  terms  of  differential  operator  equations 
(such  a;  Euler-Lagrange  equations).  Excellent  accounts  of  the  mathematical  theory  of  the  finite 
element  method  are  [Ciarlet,  1978],  [Oden  and  Reddy,  1976],  and  [Strang  and  Fix,  1973],  An 
extensive  development  from  an  engineering  point  of  view  is  presented  in  [Zienkiewicz,  1977], 

A.1 .  The  Sobolev  Spaces 

Fundamental  to  finite  element  analysis  are  a  set  of  spaces  called  the  Sobolev  spaces  (see, 
e.g.,  (Agmon,  1965],  [Adams,  1975]).  They  are  a  generalization  of  the  familiar  L2  space  which 
consists  of  all  functions  v:  fl  C  *-*  5?  (where  0  is  a  bounded  domain)  whose  L2  norm  over  fl 


is  finite.  We  denote  the  partial  derivatives  of  v  by  the  notation 


3“v(x)  = 


’The  finite  element  method  was  conceived  by  Courant  [1943). 
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where  a  —  (aIl...,a„)  is  a  multi-index  of  positive  integers  a,."  The  Sobolev  norm  of  order  m 
over  fl  combines  the  /-_>  norms  of  all  partial  derivatives  of  v  up  to  order  m: 


where  |a|  =  at  +  ■  •  •  +  an.  The  Sobolev  space  of  order  m  over  fl  is  then  defined  by 

#m(n)=  {H  IMIm.n  <  oo}.  (2) 

Clearly,  =  L2. 

We  will  also  require  the  associated  semi-norm 


which  includes  only  the  derivatives  of  order  m  exactly.  It  is  a  semi-norm  because  it  is  zero  if 
v  —  rr,u_i  €  n”l-'(n),  where  IIm— 1  is  the  space  of  polynomials  or  degree  m  —  1  defined  in  fl 
[Strang  and  Fix,  1973,  pg.  298]. 

Since  the  Sobolev  norms  are  sums  of  L-<  norms,  they  have  associated  inner  products  and, 
therefore,  the  Sobolev  spaces  are  Hilbert  spaces.3  Let  Ct?(fl)  denote  the  class  of  functions  which 
have  continuous  partial  derivatives  of  all  orders  up  to  order  q.  A  fundamental  embedding  property 
of  the  Sobolev  spaces  is  given  by  the  Sobolev  inequality  which  states  that  Cq  (Z  )t’n  )I  and  only 
if  m  —  q  >  n/2,  where  n  is  the  dimension  of  Kn. 

A. 2.  An  Abstract  Variational  Principle 

Let  V  be  a  normed  vector  space  with  norm  ||  ||,  and  S  be  a  nonempty  subset  of  V .  Moreover, 
let  a(-,  •):  V  X  V  i->  51  be  a  continuous  bilinear  form  and  / :  V  *-*  5?  be  a  continuous  linear  form. 

Definition  1  —  abstract  (quadratic)  variational  principle 
The  problem:  find  an  element  us  such  that 

us  eS  and  £(u5)  =  inf  S(vs),  (4) 

vs€S 

where  the  functional  £:V  $  is  defined  by 

f(v)  =  1q(v,  u)- /(u),  (5) 

JThe  derivatives  are  to  be  interpreted  in  the  generalized  (distributional)  sense,  but  when  a  derivative 
exists  in  the  classical  sense,  it  is  equal  to  the  generalised  derivative  (see,  e.g.,  (Schwarts,  1966]). 
Extensions  of  the  definition  of  Sobolev  spaces  and  norms  have  been  made  to  negative  and  nonintegral 
order  m  (see,  e.g.,  [Agmon,  1965],  [Adams,  1975]). 
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will  be  referred  to  as  the  abstract  variational  principle. 

Theorem  1  —  existence  and  uniqueness 
.Assume  in  addition  that 

(i)  the  space  V  is  complete,4 

(ii)  S  is  a  closed  convex  subset  of  V , 

(iii)  the  bilinear  form  a(-,  •)  is  symmetric, 

(iv)  a(-,  •)  is  V -elliptic;5  i.e.,  there  exists  a  constant  a  >  0  such  that 

Vti  (;  V,  a[v,  v)  >  o||u||2,  (6) 

then  the  abstract  variational  principle  has  a  unique  solution  us  £  S.8 

Proof.  [Ciarlet,  1978,  pg.  3]  The  bilinear  form  a(-,  •)  is  an  inner  product  over  the  space  V.  Since 
it  is  continuous,  ‘it  is  bounded,7  and  because  it  is  also  V'-elliptic,  there  exist  constants  a  and  pi 
such  that  <*|M|2  <  a(v,  v)  <  /i||u]|2.  Therefore,  the  norm  associated  with  the  inner  product  is 
equivalent8  to  the  given  norm  ||||,  and  V  becomes  a  Hilbert  space  when  it  is  equipped  with  this 
inner  product.  According  to  the  Retsz  representation  theorem  (see,  e.g.,  [Oden,  1979]  or  [Yosida, 
1971])  there  exists  an  element  u  £  V  such  that 

Vu£V,  /(u)  =  o(u,  v),  (7) 

and  because  a(-,  •)  is  symmetric, 

<f(os)  =  ^a(uS,  Vs)  —  a(u,  Vs)  =  ^a(vs  —  u,  vs  -  u)  -  ^a(u,  u). 

Therefo  e,  the  minimum  of  <f(us)  and  the  minimum  of  a(vs  —  u,  vs  —  u)  as  vs  ranges  over  the 
set  5  axe  achieved  by  the  same  element  us  £  S.  In  other  words,  solving  the  abstract  variational 
principle  is  equivalent  to  finding  an  element  us  £  S  which  iB  closest  to  u  with  respect  to  the 
norm  a(-, 

a(u  —  us,u  —  us)$=  inf  o(u  —  vs ,  u  —  vs)$.  (8) 

v*  e  s 

By  the  projection  theorem,  (see,  e.g.,  (Oden,  1979]  or  [Yosida,  1971])  the  solution  is  the  projection 
of  u  onto  S  with  respect  to  the  inner  product  o(-,  •),  and  its  existence  and  uniqueness  is  assured 
by  the  fact  that  S  is  a  closed  convex  subset  of  V.  | 

4 That  is  to  say,  it  is  a  Banach  space. 

5  V-ellipticity  means  that  the  bilinear  form  is  positive  definite;  i.e.,  o(ti,  v)  =  0  if  and  only  if  v  —  0. 
6Theorom  1  is  a  generalisation  of  the  familiar  theorem  for  the  existence  of  a  unique  solution  to  a 
(quadratic)  minimization  problem  in  mathematical  programming  (see  e.g.  |Lucnberger,  1973]). 

7 A  bilinear  form  is  continuous  if  and  only  if  it  is  bounded;  i.e.,  there  exists  a  constant  pi  such  that 
|o(u,  v)|  <  mIMIIMI  (Rektorys,  1980,  pg.  111]. 

"Two  norms  ||-||  and  |H|'  on  a  linear  vector  space  V  are  called  cquivsle.t  if  the  corresponding 
metrics  are  equivalent.  This  amounts  to  the  existence  of  two  positive  constants  c,  and  ci  such  that 

ciiN  <  IN'  <  «*IM|. 
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The  abstract  formulation  encompasses  linear  variational  problems  which  are  posed  classically 
in  terms  of  variational  principles  involving  the  minimisation  of  quadratic  functionals  £ (v)  = 
^a(v,  v)  —  /( v)  over  an  admissible  space  of  Junctions  v  6  V.  Such  functionals  often  represent  the 
potential  energy  of  a  physical  system,  a(v,  v)  being  the  second-degree  term  which  is  the  strain 
energy  in  the  function  v  ( /(v )  is  a  first-degree  term).  The  associated  inner  product  a(v,  to)  is  the 
energy  inner  product  which  is  intrinsic  to  the  particular  variational  principle,  and  is  defined  for 
all  admissible  functions  v  and  to. 

It  is  clear  from  the  above  discussion  that  the  admissible  space  V'  must  be  complete  and  that 
£[v)  must  be  well-defined  for  all  v  £  V  (i.e.,  that  V  must  be  a  space  of  finite  energy).  The  Sobolev 
spaces  fulfill  these  conditions.  Their  use  as  generalized  energy  spaces  is  natural  in  the  sense  that, 
for  a  given  variational  principle,  the  energy  inner  product  iB  well-defined  over  a  Sobolev  space 
lfm,  where  m  is  the  highest  order  of  partial  derivative  of  v  which  occurs  in  a(v,  v).  In  genera) 
then,  V  is  the  space  M m  whose  natural  norm  is  the  Sobolev  norm  ||-||m. 

The  role  of  the  subset  S  is  in  approximating  the  exact  solution  u.  Although  u  is  usually 
impossible  to  obtain  over  the  full  admissible  space  V,  it  may  be  relatively  straightforward  to 
optimally  approximate  it  by  an  element  us  £  S,  especially  if  5  is  taken  to  be  a  closed  subspace 
cf  V.  The  approximation  is  optimal  in  the  sense  of  equation  (8).  In  the  ensuing  discuBBion,  we 
will  restrict  ourselves  to  the  special  case  where  S  is  a  closed  subspace  of  V.  The  approximate 
solution  of  the  abstract  variational  principle  may  then  be  characterized  by  the  following  theorem. 

Theorem  2  —  variational  equation 

If  S  is  a  closed  subspace  of  V,  then  us  £  S  is  o  solution  of  the  abstract  variational  principle 
if  and  only  if  it  satisfies  the  variational  equation 

Vvs  G  S,  «(us,  vs)  =  /(vs).  (9) 

Proof.  If  tts  minimizes  £  over  S,  then  for  any  e  and  vs  6  S, 

£ (tts)  <  £(us  +  £vs)  =  ^o(us  +  tvs ,  u5  -f  evs)  —  f[ us  +  «vs) 

=  £(us)  +  <[<*(«*,  Vs)  —  /(vs)]  -I-  ^«Sa(vS,  Vs). 

Therefore, 

0  <  e[a(us,  vs)  —  /(vs)]  +  ^e2a(vS>  vS), 

and  since  this  must  be  true  for  small  e,  both  positive  and  negative,  it  follows  that  o(us,  Vs)  = 

/(vs).9  I 

*In  the  general  case  where  S  is  not  a  closed  subspace,  but  is  only  a  closed  convex  subset  of  V 
as  required  by  condition  (ii)  of  Theorem  1,  the  solution  us  must  satisfy  the  vsrtahonoJ  inequality 
a(us,  vs  -  us)  >  /(vs  -  vs)  (Ciarlet,  1978,  pg.  3). 
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We  will  now  discuss  some  important  properties  of  the  solution  of  the  abstract  variational 
principle.  First,  from  the  proof  of  the  theorem,  it  is  clear  that  (9)  is  the  well-known  condition 
for  the  vanishing  of  the  first  variation  of  <  at  us,  in  the  direction  of  vs.  In  particular,  if  S  is  the 
whole  space  V,  then  the  solution  satisfies 

o(«.  «)  =  f{v),  (10) 

and  the  first  variation  at  u  vanishes  in  every  direction  v.  Setting  u  =  v,  we  have  a(u,  u)  =  /( u) 
and  hence 

f(u)  =  ^o(u.  u)  -  /(«)  =  —  «); 

i.e.,  o£  the  minimum,  the  strain  energy  is  the  negative  of  the  potential  energy.  Now,  E (u)  <  £ (us) 
since  u  is  minimal  over  a  wider  class  of  functions.  Then, 

a(us,  u5)  <  a(u,  u), 

and  so  the  strain  energy  in  us  always  underestimates  the  strain  energy  in  u.  Moreover,  since  us 
is  the  projection  of  u  onto  the  subspace  S,  the  error  es  =  u  —  us  is  orthogonal  to  S: 

a(es,  us)  =  0  Vvs€  S. 

In  particular,  o(es,  us)  —  0  or  a(ti,  us)  —  a(tzs,  us)  and 

o(es,  es)  =  a(u,  u)  —  o(ns,  us ); 

i.e.,  the  energy  in  the  error  equals  the  error  »n  the  energy  (the  Pythagorean  theorem  holds). 

A. 3.  The  Ritz  Approximation 

The  key  condition  in  the  hypothesis  of  Theorem  2  is  that  the  subspace  S  be  closed.  How  can 
we  ensure  this?  One  possibility  arises  from  the  fact  that  finite  dimensional  subspaces  are  always 
closed.  A  number  of  classical  methods  for  solving  variational  problems,  called  direct  methods,  are 
based  on  this.10  One  of  these,  the  (Raleigh-)  Ritz  method  [Ritz,  1908;  Mikhlin,  1964;  Rektorys, 
1980]  is  of  fundamental  importance  when  a  variational  principle  is  involved.  In  the  Ritz  method, 
we  chocae  a  finite  dimensional  subspace 

N 

s  =  sh  ~  {vh  1  v*  =  v><>  (H) 

i— i 

where  p, . are  independent  basis  functions  which  span  Sh  and  Vi . v/v  are  unknown 

real  parameters. 

10 Direc*  methods  include  the  metho.J  of  weighted  residuals  whose  special  cases  include  collocation  methods 
and  the  (jaierkin  method ,  the  method  of  orthonormal  (t.g.,  Fourier)  series ,  and  the  least  squares  method,  (see, 
e.g.,  [FinUyaon,  1972),  [Mikhlin,  1964),  [Rektorys,  1980),  [Zienkiewics,  1977)). 
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Au  =  f,  (13) 

where  A  G  —  [a(<£,,  <f>j)\  and  where  f  £  SUN  =  [/(<£,)]  and  u  £  KN  =  |u,],  called  the  discrete 

variational  equation.  Since  the  matrix  of  coefficients  A  is  nonsingular,  the  discrete  solution  is 
given  by  u  —  A- 1  f ,  although  for  the  problem  at  hand  A  is  huge,  so  it  is  usually  unpractical 
to  compute  A~'  directly.  In  the  next  section,  we  describe  special  types  of  basis  functions  which 
ultimately  lead  to  practical  iterative  solutions. 


A. 4.  Finite  Element  Spaces 

The  Ritz  method  has  given  us  a  discrete  solution  uh  =  u,^  which  is  optimal  in  the 

sense  that  the  energy  in  the  error,  as  measured  in  the  natural  energy  norm  a(u  —  uh,  u  —  uA)i, 
is  as  small  as  possible.  In  the  classical  Rita  method,  the  basis  functions  4>,  are  generally  chosen 
to  be  fairly  complicated  functions  which  have  global  support  over  the  domain  in  question  (e.g., 
trigonometric  functions)  [Mikhlin  1964;  Rektorys,  1980].  Although  this  choice  may  be  beneficial 
for  analytic  purposes,  it  renders  the  method  unsuitable  for  numerical  computation.  The  problem 
is  overcome  by  the  finite  element  method  which  is  a  systematic  procedure  for  constructing  finite 
dimensional  approximating  subspaces  Sh,  called  finite  element  spaces,  which  are  very  convenient 
for  numerical  computation.  In  certain  forms,  the  method  may  be  considered  to  be  a  special 
instance  of  the  Rits  method  in  which  the  basis  functions  are  simple  functions  having  local  support. 
In  the  ensuing  discussion,  we  will  restrict  ourselves  to  a  domain  fl  which  is  a  polygon  in  Ut2  with 
boundary  dfl  .**  The  following  are  basic  characteristics  of  the  construction  in  its  simplest  form: 

(i)  A  * triangulation *  Th  is  established  over  the  domain:  H  =  UesT*  that  is,  the  domain 
is  partitioned  into  the  union  of  subdomains  E  E  Th  called  finite  elements,  such  that  the  E 
are  closed  sets  with  nonempty  interiors  and  polygonal  boundaries.  TLj  elements  are  usually 

"The  theory  has  been  extended  to  domains  with  curved  boundaries  in  any  number  of  dimensions. 
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adjacently-placed  triangles  or  rectangles  which  overlap  only  at  the  inter-element  boundaries. 
Associated  with  the  triangulation  is  its  fundamental  length  h.12 

(ii)  The  elements  are  considered  to  be  interconnected  at  a  discrete  number  of  points  on  the 
inter-eleinent  boundaries  which  are  called  the  nodes  of  the  triangulation.  The  unknown  real 
parameters  of  the  discrete  problem  are  the  nodal  variables,  the  values  of  the  solution  (and/or 
possibly  of  its  derivatives)  at  the  nodes. 

(iii)  Associated  with  the  tri angulation  is  a  space  of  functions  5h  defined  over  Cl.  Defined  within 

each  element  E  is  a  finite  dimensional  space  PK  =  (vh|K  |  €  Sh}  consisting  of  functions 

which  are  polynomials  (or  ratios  of  polynomials).  The  polynomials  pE  £  PK  represent  a  local 
approximation  of  the  solution  within  E,  and  arc  uniquely  determined  by  the  nodal  variables 
associated  with  the  element. 

(iv)  In  certain  elements,  the  functions  pE  may  have  to  satisfy  the  essential  boundary  conditions 
of  the  problem. 

While  the  classical  Ritz  method  is  limited  to  geometrically  simple  domains  Cl,  in  the  finite 
element  method  this  limitation  occurs  only  within  the  element  itself.  Consequently,  it  is  possible 
to  “assemble”  complicated  configurations  from  simple  element  shapes.  Several  factors  contribute 
to  the  success  of  the  finite  element  method  from  a  computational  point  of  view.  Firstly,  due  to 
the  fact  that  in  the  element  interiors  the  solution  is  approximated  by  a  low-order  polynomial  in 
x  and  y,  the  computations  required  to  compute  the  discrete  functional  in  (12)  or,  equivalently, 
to  compute  the  entries  of  matrix  A  of  the  discrete  variational  equation  (13)  are  often  Biinpie. 
Secondly,  it  can  be  shown  that,  associated  with  the  local  polynomial  functions,  there  exists  a 
canonical  set  of  basis  functions  <(>,  spanning  Sh  which  sire  also  piecewise  polynomials  and  which 
have  local  support.  A  will  therefore  be  sparse  and  banded;  that  is,  most  of  its  relatively  few 
nonzero  entries  will  he  near  the  main  diagonal.  Thirdly,  when  the  problem  is  well-posed  in  terms 
of  a  variational  principle,  a(-,  •)  will  be  symmetric  and  ^''-elliptic,  which  guarantees  that  A  will 
be  nonsingular,  symmetric,  and  positive  definite.  In  addition  to  these  clear  merits,  piecewise 
polynomials  are  remarkable  in  that  they  are  optimal  in  terms  of  their  approximation  properties 
and  in  that  these  properties  are  essential  for  proving  convergence  of  the  method  [Strang  and  Fix, 
1973,  pg.  153]. 

The  convergence  properties  of  the  finite  element  method  are  an  important  issue.  The  object 
of  the  Ritz  method  is  to  find  optimal  values  for  the  nodal  variables  (which  are  the  parameters  of 
the  discrete  solution)  by  minimizing  the  discrete  functional  £(vE).  This  suggests  immediately  the 
possibility  of  approximating  the  exact  solution  u  by  a  minimising  sequence  of  discrete  solutions  to 
discrete  problems  associated  with  a  family  of  subspaces  Sh  whose  fundamental  length  h  has  limit 
zero.  Although  the  approximation  is  known  by  (8)  to  be  optimal  in  terms  of  the  norm  a(-,  -)i,  it 
is  more  convenient  to  analyze  the  error  in  terms  of  the  natural  Sobolev  norm  ||-|jm  of  V  £  l/m. 
The  following  theorem  gives  a  sufficient  condition  for  the  convergence  of  such  a  sequence. 

12The  fundamental  length  h  of  the  triangulation  Th  is  the  maximum  “radius”  of  the  elements.  As  the 
subdivision  is  made  finer,  the  number  of  elements  increase  and  h  — »  0. 
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Theorem  3  —  Cea’s  Lemma 

Since  there  exists  a  constant  C  independent  upon  the  subspace  Sh  such  that 


then  a  sufficient  condition  for  convergence  is  that  there  exists  a  family  SK  of  subspaces  of  the 
space  V  such  that,  for  each  u  £  V ,  lim/l_o||u  —  u,‘||  =  0;  i.e., 

litn  inf  ||u  —  v/lll  =  0. 
fc— o  v'es* 


Proof.  Equation  (14)  follows  from  (8)  due  to  the  continuity  and  V-ellipticity  of  o(-,  •).  Moreover, 
C  —  \J\  is  a  constant  independent  upon  the  subspacc  Sh.  | 

We  see  then  that  an  estimation  of  the  error  reduces  to  finding  the  distance  between  the  exact 
solution  u  and  the  subspace  Sh  —  a  problem  in  approximation  theory.  The  basic  hypothesis  about 
the  finite  element  space  Sh  was  that  the  finite-dimensional  space  PE  within  each  element  E  is  a 
space  of  polynomials.  If  we  assume  that  the  space  contains  the  complete  polynomials  of  degree  k 
(i.e.,  II*  C  P1'),  it  can  be  shown  in  general  that  the  approximation  error  in  the  s1*1  derivative, 
where  m,  is  of  the  form 


||u  -  uh||,  =  0(hfc+I~s  +  h2(fc+1“m))  (15) 

[Strang  and  Fix,  1973,  pg.  107].  On  the  other  hand,  because  the  approximation  minimises  the 
strain  energy,  the  order  of  convergence  of  the  derivative  is  better.  It  is  order  m). 

The  convergence  properties  implied  by  Cla’s  lemma  are  contingent  upon  the  finite  element 
spaces  Sh  being  subspaces  of  the  admissible  space  V.  In  view  of  this,  if  the  energy  inner  product 
a(v.  v)  involves  partial  derivatives  of  v  of  order  m  so  that  V  C  J/m(fi),  ensuring  convergence 
amounts  to  imposing  the  following  two  requirements  on  the  local  functions  pE: 

(i)  Completeness  condition:  As  the  size  of  any  element  tends  to  zero,  the  function  pE  must  be 
such  that  a  constant  value  of  the  m'*1  derivative  will  be  attainable  in  the  element  subdomain; 
i.e.  we  must  have  k  >  m  so  that  PE  C  Mm(E),  VE  G  Th. 

(ii)  Conformity  condition:  All  derivatives  of  pE  of  order  less  than  m  must  be  continuous  across 
inter-element  boundaries;  that  is,  Sh  C  Cm~ *(fl). 

The  two  requirements  are  necessary  and  sufficient  for  Sh  C  l/m(fi)  when  the  pE  are  polynomials 
(or  ratios  of  polynomials).  Another  way  of  stating  the  completeness  condition  is  that  the  local 
polynomials  must  be  able  to  reproduce  a  state  of  constant  strain  —  any  solution  which  is  a 
polynomial  of  degree  m.  When  the  local  polynomials  satisfy  the  conformity  condition,  the  elements 
are  called  conforming  finite  elements,  and  their  use  leads  what  are  referred  to  as  conforming  finite 
element  methods. 
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A.S.  Nonconforming  Elements 

In  the  above  discussion,  we  assumed  that  conforming  finite  element  methods  approximate 
the  solution  u  of  £(u)  =  inf^gs  £(v)  by  the  solution  uh  of  £(uh)  —  infvke$>  <f  (vh)  where  S'* 
is  a  subspace  of  S.  This  is  a  global  condition  on  the  approximation  which  is  often  violated  for 
reasons  of  computational  convenience.  For  instance,  it  may  be  violated  by  dropping  the  element 
conformity  condition.13  Elements  which  do  so  are  called  nonconforming  elements.  They  are  often 
used  in  practice  for  higher-order  problems  because  conforming  elements  tor  such  problems  are 
unnecessarily  complicated  or  must  have  a  large  number  of  degrees  of  freedom  in  order  to  satisfy 
the  inter-element  conformity  conditions. 

If  nonconforming  elements  are  used,  it  is  clearly  impossible  to  evaluate  the  true  energy 
functional  £(vh),  due  to  the  singularities  in  the  mtc  derivatives  of  vh  which  occur  at  the  element 
boundaries.  To  avoid  this  problem,  we  can  simply  ignore  the  discontinuities  between  elements  by 
computing  the  strain  energies  within  each  element  and  then  summing  the  individual  contributions; 
that  is,  for  the  original  energy  inner  product  a(-,  -),  we  substitute  the  approximate  energy  inner 
product  aA(-,  ■)  defined  by  the  bilinear  form 

«*(•.  •)  =  E  4.  -)\b,  (»6) 

A'€T» 

where  the  notation  jt;  means  a  restriction  to  the  element  domain.  The  approximate  variational 
principle  is  then  the  problem  of  finding  a  u h  £  Sh  which  minimises  the  functional 

&(**)-  •*)-/(•*).  (17) 

and  the  necessary  condition  for  the  vanishing  of  the  first  variation  becomes 

Vvh  6  S\  ah(u\  vh)  =  /(o'*).  (18) 

Following  in  the  spirit  of  the  conforming  case,  we  must  determine  sufficient  conditions  for 
the  existence  and  uniqueness  of  the  solution  uh  to  the  approximate  variational  principle,  as  well 
as  under  what  conditions  this  approximate  solution  converges  to  the  exact  solution  u  as  h  -»  0. 
The  conditions  are  natural  extensions  of  Theorem  1  and  Cda’s  lemma,  and  are  given  in  the 
following  theorem. 

Theorem  4  —  existence  and  uniqueness  (nonconforming  case) 

If 

(i)  there  exists  a  mapping  ||-|j(,:S'1  •-»  R  which  is  a  norm  over  S'*, 

(ii)  Oh(-,  •)  is  bounded  and  Sh-elliptic,  in  that  there  exists  a  constant  an  >  0  such  that 

13This  violation  is  an  example  of  a  so  called  variations!  crime  [Strang  and  Fix,  1973,  Chapter  4). 
Besides  violation  of  the  element  conformity  condition,  variational  crimes  also  inc'ude  inexact  evaluation 
of  the  functional  £(v*)  (i.e.,  of  the  quadratic  form  a(u\  v*)  and  linear  form  /(v*1))  such  as  by  numerical 
integration,  as  well  as  various  techniques  for  the  approximation  of  essential  boundary  conditions. 
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6  s’1,  atl(vA,  t/*)  >  afc||«*|(J, 

then  the  approximate  variational  principle  has  a  unique  solution  uh  G  SK. 

Proof.  Refer  to  the  discussion  in  [Ciarlet,  1978,  Chapter  4].  | 

On  the  other  hand,  to  obtain  convergence,  we  impose  a  stronger  condition,  uniform  S  Ihpticity, 
which  requires  that  there  exist  a  constant  a  >  0,  independent  of  h,  such  that 


VS*,  Vv*  e  S\  ah{vh,  vh)  >  a||«fc||J. 
Convergence  is  then  guaranteed  by  the  following  theorem. 


(19) 


Theorem  5  —  Strang’s  lemma 

Given  a  family  of  discrete  problems  for  which  the  associated  approximate  energy  inner 
products  are  uniformly  Sh-elliptic,  then  there  exists  a  constant  C ,  independent  of  S’1,  such  that 


11“  -  “hIU  <  Cl  inf  ||u  -  v*||fc  +  si 

i wk 


sup 

6Sk 


jaA(u,  wh)  —  /( 


l|whlk 


•^V 


(20) 


Proof.  See  [Ciarlet,  1978,  pg.  210].  | 


Strang’s  lemma  is  a  generalization  of  Cea’s  lemma  for  conforming  elements  —  in  addition 
to  the  usual  approximation  error  term,  we  have  the  (inf)  term  which  measures  the  consistency 
error  of  the  nonconforming  space.  Since  the  difference  Oh(u,  wh)  —  f{u>h)  is  zero  for  all  wh  E  Sh 
when  Sh  C  V,  the  consistency  error  for  conforming  spaces  is  identically  zero  and  Strang’s  lemma 
reduces  to  Cea's  lemma.  However,  for  the  nonconforming  case,  convergence  is  obtained  if  the 
consistency  condition, 


lim  sup 

k—0 


K(t*.  uh)  —  /(w'*)| 

IKIk 


=  o, 


Vv*  6  Sh, 


(21) 


is  satisfied. 


The  consistency  condition  was  first  recognized  empirically  and  was  stated  in  the  form  of  a 
simple  test  known  as  the  patch  test.  Subsequently,  Strang  proved  mathematically  that  the  patch 
test  was  indeed  a  test  for  consistency  and,  by  essentially  making  the  above  arguments,  that 
nonconforming  elements  which  pass  it  will  yield  convergence  (see  [Strang  and  Fix,  1973,  Chapter 
4]).  The  test  remains  a  convenient  one  to  apply. 

Theorem  6  —  the  patch  test 

Suppose  that  the  energy  inner  product  a>,(u,  v)  contains  derivatives  of  order  at  most  m  and 
the  nonconforming  space  Sh  contains  all  polynomials  irm  of  degree  at  most  m.  If  the  nonconforming 
finite  element  method  recovers  all  solutions  which  are  polynomials  of  degree  at  most  m,  then  the 
patch  test  is  passed,  and  lim^— >ollu  —  —  0. 
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In  other  words,  suppose  that  we  put  an  arbitrary  patch  of  nonconforming  elements  associated 
with  the  nonconforming  space  Sh  in  a  state  of  constant  strain;  that  is,  we  impose  vK  =  wm  £  nm 
on  the  displacements  at  nodes  around  the  patch  boundary.  Because  the  completeness  condition 
of  the  previous  section  is  still  binding  on  nonconforming  elements,  this  polynomial  is  both  an 
element  of  Sh  and  an  element  of  V,  hence  its  consistency  error  is  sero.  The  conforming  (Riti) 
solution  to  (12)  or  (13)  and  the  nonconforming,  discrete  solution  of  the  approximate  variational 
principle  ought  then  to  be  identical  and  equal  to  rm.  The  test  is  to  determine  whether  this  is 
indeed  the  case. 
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B.  ITERATIVE  SOLUTION  OF  LARGE  LINEAR  SYSTEMS 

The  approximation  of  a  variational  principle  by  direct  methods  such  as  the  finite  clement 
method  (or  the  approximation  of  a  boundary  value  problem  by  finite  differences)  gives  rise  to  a 
system  of  simultaneous  algebraic  equations.  For  quadratic  functionals  (or  linear  boundary  value 
problems),  the  system  will  be  linear.  In  this  chapter,  we  consider  the  problem  of  solving,  by 
iterative  means,  systems  of  linear  equations  of  the  form 

N 

a')u> =  /■>  (i) 

j- 1 

where  the  coefficients  a,}  and  the  va'ues  /,  are  known.  The  system  may  also  be  written  as  the 
matrix  equation  _ 

Au  =  f,  (2) 

where  A  £  fRNN  =  [a,j]  is  a  nonsingular  matrix  and  where  f  6  =  [/,]  is  a  column  vector.  We 

wish  to  solve  the  system  for  the  column  vector  u  £  fkN  =  [u,]  =  A~'f.  In  applying  the  finite 
clement  method  to  our  visual  surface  reconstruction  problem,  we  will  obtain  large  sparse  matrices 

A.  In  this  appendix  we  will  be  concerned  with  numerical  methods  which  are  capable  of  solving 
such  systems  where  N  is  in  the  range  of,  say,  103  —  106,  or  larger. 

An  iterative  method  for  solving  equations  (2)  computes  a  sequence  of  approximations 
u(i),  u(2), . . . ,  to  the  exact  solution  u.  A  new,  and  hopefully  better,  approximation  is  computed 
at  each  iteration,  but,  in  general,  the  exact  solution  cannot  be  obtained  in  a  finite  number  of 
iterations.  If,  regardless  of  the  initial  approximation  u*°V4  the  approximation  error  (i.e.,  the 
difference  between  the  exact  solution  and  the  approximations,  measured  in  some  appropriate  norm) 
tends  to  zero  as  the  number  of  iterations  increases,  the  iterative  method  is  said  to  be  convergent, 
and  the  rate  at  which  the  approximation  error  tends  to  zero  is  called  its  rate  of  convergence.  In 
order  that  an  iterative  method  be  of  practical  use,  it  is  important  that  it  be  convergent,  and  that 
it  exhibit  a  sufficiently  large  rate  of  convergence  to  an  approximation  of  prespecified  accuracy. 
In  this  appendix,  we  review  a  number  of  the  most  common  iterative  methods,  and  examine  their 
convergence  properties  and  their  rates  of  convergence.  References  for  this  material  are  [Forsythe 
and  Wasow,  I960],  [Smith,  1977],  [Varga,  1963]  and  [Young,  1971;  Young  and  Gregory,  1972, 
1973]. 

B. l .  Basic  Relaxation  Methods 

Let  us  assume  that  A  has  nonzero  diagonal  elements.  It  will  be  convenient  to  express  A  as 
the  sum 

A  =  D  —  L  —  U, 

MThe  trivial  initial  approximation  u(0)  =  0  is  usually  chosen. 
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where  D  £  ?RNN  is  diagonal,  L  £  $tNN  is  strictly  lower  triangular,  and  U  €  RNN  is  strictly  upper 
triangular.  Clearly  the  equations  in  (1)  can  be  rewritten  in  the  following  form: 

<*»«,  —  —  5^  oi;u,  +  /„  1  <  *  <  N. 

i  <jSN 

rV* 

Next,  we  will  define  three  basic  iterative  methods  for  solving  (2),  popularly  known  as  relaxation 
method s  in  the  context  of  the  numerical  solution  of  partial  differential  equations. 

Jacobi  Relaxation 

The  Jacobi  method  (or  the  method  of  simultaneous  displacements)  is  defined  by 

4*  fit  1  <  i <  N, 

which  can  be  written  in  matrix  form  as 

Du(fc+"  =  (L  +  U)u,lc*  +  f, 

thus  giving  us  the  iterative  scheme 

„(*+i)  =  D-1(L4-U)u<fc>4~D-,f.  (3) 

The  matrix 

Gj  =  D~*(L  4-  V) 

is  called  the  Jacobi  itc.ation  matrix  associated  with  matrix  A. 

Clearly,  the  Jacobi  method  is  a  parallel  method  because  the  elements  of  the  new  approximation 
u(fc+1)  may  be  computed  simultaneously  and  in  parallel  by  a  network  of  processors  whose  inputs 
are  elements  of  the  old  approximation  u^.  As  such,  it  requires  the  storage  of  both  the  old  and 
the  new  approximations. 

Gauss-Seidel  Relaxation 

Convergence  of  the  Jacobi  method  is  usually  very  slow.  In  a  closely  related  method,  the 
so  called  Gauss-Seidel  method  (or  the  method  of  immediate  displacements),  the  elements  of  the 
new  approximation  are  used  in  subsequent  computations  immediately  after  they  become  available. 
This  increases  the  rate  of  convergence  somewhat,  but  typically  less  than  an  order  of  magnitude. 
The  equations  are  written  as 

_  -  '£;  „«}»+■'  -  £  „„«!*> +/.,  i  < ,  <  w, 
i-i  r-H-i 

which,  in  matrix  form,  becomes 

Du(fe+i)  =  Lu<*+,)  4-  Uu(fc)  4-  f,  (4) 
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from  which  we  obtain  the  iterative  scheme  defined  by 

„<*+')  =  (D-L)-1Uu(*r)-t-(D-L)_1f.  (5) 

The  Gauss-Seidel  iteration  matrix  associated  with  matrix  A  is  therefore  given  by 

Gas  =  (D-L)-'U. 

In  the  Gauss-Seidel  method,  we  no  longer  need  to  store  simultaneously  the  old  and  new 
approximations.  As  they  arc  computed,  the  new  elements  can  simply  displace  the  old  ones. 
Moreover,  since  the  new  values  are  exploited  immediately  in  subsequent  computations  we  can 
intuitively  expect  a  higher  rate  of  convergence  compared  with  the  Jacobi  method.  On  the  other 
hand,  it  can  easily  be  seen  that  the  Gauss-Seidel  method  is  inherently  a  sequential  method  which 
renders  it  unsuitable  for  implementation  as  a  parallel  network. 

The  Successive  Overrelaxation  (SOR)  Method 

The  convergence  of  the  Gauss-Seidel  method  may  be  accelerated  by  a  simple  modification. 
Let  us  define  the  dynamic  residual  vector 15  at  the  k ^  iteration  of  the  Gauss-Seidel  relaxation 
method  as 

r(fc)  _  „(*+!)_„<*) 

=  D-,(Lu<*+,)  4  Uu,k)  +  f)  -  u<k> 

(see  (4)).  Then,  it  is  evident  that  the  Gauss-Seidel  scheme  can  be  written  in  the  form 

„(*+i)  =  u(fc)  +  r(*). 

If  the  elements  of  all  successive  residual  vectors  are  one-signed  (as  they  usually  are  when 
approximating  elliptic  problems),  then  it  is  reasonable  to  anticipate  an  acceleration  in  the 
convergence  of  the  Gauss-Seidel  scheme  if  the  residual  vector  is  scaled  by  a  fixed  real  number 
w  >  1  before  it  is  added  to  u^  in  each  Gauss-Seidel  iteration.  u>  is  called  the  relaxation  parameter. 
This  idea  leads  to  one  of  the  most  successful  iterative  methods  for  solving  systems  of  linear 
equations  currently  available,  the  successive  overrelaxation  method.  The  method  may  be  defined 
by 

u(*+i)  = 

=  u<*>  4  w[D-,(Lu(k+1>  -(-  Uu(k)  -f  f)  -  u(k)], 
which  we  can  manipulate  into  the  form 

u(fc+>)  =  (I  -  wD— lL)~ 1  ({I  -  w)I  +  wD-‘U]u<k>  4  (I  -  wD-,L)~1wD-,f.  (6) 

The  successive  overrelaxation  iteration  matrix  is  therefore  given  by 

Gw  =  (I  —  wD— 'L)— '[(1  —  o/)I4  wD— 'U). 

,sThe  residual  vector  is  also  called  the  correction  or  displacement  vector.  Note  that  the  residual  of  an 
equation  is  the  amount  by  which  the  equation  fails  to  be  satisfied. 
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Three  cases  arise: 

(i)  if  u  >  1  the  scheme  is  called  overrelaxation,16 

(ii)  if  w  =  1  we  obtain  the  Gauss-Seidel  method,  and 

(iii)  if  0  <  w  <  1  the  scheme  is  termed  underre taxation. 

Conditions  for  Convergence 

The  Jacobi,  Gauss-Seidel,  and  successive  ovcrrelaxation  methods  are  particular  instances  of 
the  general  stationary  iterative  method 

u^UGuW+e,  (7) 

where  the  iteration  matrix  G  is  taken  to  be  G j,  Gas,  and  Gw  respectively,  and  c  is  a  known 
column  vector.17  By  subtracting  u  =  Gu  +  e  from  (7),  wc  can  obtain  an  expression  for  the  errors 
=  u'*>  —  u  of  successive  approximations, 

«(*+»  =  Ge!t> 

=  G*+,e(0>.  (8) 

The  sequence  of  approximations  converges  to  the  solution  u  if  liroit—oo  =  0,  which  will 

clearly  be  the  case  if  and  only  if  lim/f-^ooG*  =  0,  since  u,0l  (and  hence  e(0))  is  arbitrary.  Let 
G  have  eigenvalues  Xi ,  X2,  • . . ,  X^,  and  assume  that  the  corresponding  eigenvectors  v, ,  Va, . . .  ,v,v 
are  linearly  independent.  Then  the  initial  error  vector  can  be  expressed  uniquely  as  the  linear 
combination 

N 

C(0)  ^  C,  V  j  • 

1  =  1 

But,  by  (8), 

e<*>  =  GktW 
N 

~  51  e'Gky< 

.  =  1 

N 

—  53  c,\kvt. 

•  «  1 

It  follows  that  in  the  limit,  e*fc*  will  be  zero  for  an  arbitrary  initial  approximation  if  and  only  if 
|X,|  <  1,  for  1  <  »  <  N.  Thus,  we  have  the  following  theorem. 

"’ll  should  be  noted  that  there  exist  classes  of  matrices  (which  arise  from  many  first  and  second 
order  partial  differential  equations)  for  which  the  optimal  value  of  u>,  yielding  the  largest  rate  of 
convergence,  may  be  determined  analytically  (see  e.g.  (Young,  1972),  (Young  and  Gregory,  1973]).  Often, 
the  convergence  may  be  adequately  accelerated  by  not  necessarily  optimal  values  of  u  chosen  empirically. 
Of  course,  it  is  possible  to  vary  u  from  iteration  to  iteration  or  from  one  equation  to  the  next.  A 
number  of  these  modified  methods  have  been  studied  in  the  literature  (see  e.g.  [Varga,  1962]  or  [Young, 
1971;  Young  and  Gregory,  1972,  1973[). 

17  The  method  is  obtained  by  writing  the  original  system  Au  =  f  as  a  =  Gu  +  e,  and  is  referred  to  as 
being  stationary  because  G  is  fixed  for  all  iterations. 
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Theorem  1  - —  necessary  and  sufficient  condition  for  convergence  of  the  stationary  iterative  method 

The  stationary  iterative  method  (7)  is  convergent  if  and  only  if 

p(C.)  =  niax|X,(G)|  <  1, 

i=l 

where  p( G)  is  called  the  spectral  radius  of  G.19 

Theorem  1  is  mainly  of  theoretical  value  because,  in  practice,  it  is  difficult  to  determine 
the  eigenvalues  of  G.  Fortunately,  we  have  a  useful  corollary  giving  sufficient  conditions  for 
convergence.  The  corollary  results  from  the  observation  that  for  some  matrix  norm  ||G||  which  is 
consistent19  with  a  vector  norm  | j v 1 1 j , 

l|e'fc>il  <  l|Gfc||  ||e‘°>||  <  ||G||fc||e<°)||.  (9) 

Hence  we  obtain  the  following  corollary  to  Theorem  1. 


Corollary  1  —  sufficient  condition  for  convergence  of  the  stationary  iterative  method 
//  1 1 G 1 1  <  1,  then  the  stationary  iterative  method  (7)  ts  convergent.20 

Suppose  that  the  stationary  iterative  method  is  convergent.  It  can  be  shown  (see  [Varga,  1962]  or 
[Young,  1971])  that 

Km  \\Gk\\,.2'/k  =  P(G). 

fc— oo 

Hence,  from  (9)  we  have  that,  for  large  k, 

\\''k)\\L,  ~  p{G)k\\'W\\Lt. 

Thus,  in  a  certain  sense,  p(G)  is  a  measure  of  the  rate  of  convergence  of  the  iterative  method  and 
therefore,  like  convergence  itself,  depends  on  the  eigenvalues  of  G. 


We  illustrate  an  application  of  Corollary  1  in  obtaining  a  simple  but  important  sufficient 
condition  for  the  convergence  of  Jacobi  relaxation.  The  Jacobi  iteration  matrix  G  j  consists  of 
the  elements 


9u  ~  0. 


Therefore  the  L^,  norm  of  Gj  is  given  by 


l!Gj||oo  =  max  53 


,8It  can  be  shown  (see  e.g.  [Varga,  1962])  that  the  theorem  holds  without  the  independence  assumption 
on  the  eigenvectors. 

'®If  a  matrix  norm  and  a  vector  norm  satisfy  the  relation  ||Au||  <  ||A||  ||u||  for  any  A  and  u,  then  the 
two  norms  are  said  to  be  consistent  [Dahlquist  and  Bjorck,  1974,  pg.  175]. 

20The  condition  is  not  a  necessary  one  because  we  can  have  the  case  that  ||G||  >  1  when  p(G)  <  1. 
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Clearly,  if  |a„|  >  |a«j|.  ||Gj||oo  <  1,  and  Jacobi  relaxation  will  converge  by  corollary 

1.  A  much  more  general  result  may  be  obtained.  We  begin  by  defining  two  important  properties 
that  A  may  possess. 

Definition  1  —  (weakly)  diagonally  dominant  matrix 

A  matrix  A  of  order  N  is  said  to  be  diagonally  dominant  if 

i°tti  >  E  K>i.  !<•<*■ 

1  SjSW 
)¥=' 

The  matrix  is  said  to  be  weakly  diagonally  dominant  if  the  >  relation  in  the  above  inequality  can 
be  replaced  by  >  m  some,  but  not  all,  oj  the  equations. 

Definition  2  —  irreducible  matrix 

A  matrix  A  of  order  N  is  irreducible  if  and  only  if  N  =  1,  or  if  N  >1  and  for  any  i 
and  j  such  that  1  <  i,j  <  N  and  i  =£  j  either  atJ  =  0,  or  there  exist  ij,  s'j, . . . ,  i*  such  that 

a,,,,  a,  1)t,  i .  .a,k tj  0. 

It  can  be  shown  that  if  A  is  irreducible  and  has  weak  diagonal  dominance,  then  it  is  nonsingular, 
and  if  in  addition  it  is  symmetric  and  has  non-negative  diagonal  elements,  then  it  is  positive 
definite.  The  general  theorem  is  given  next  (For  a  proof  see  [Varga,  1962,  pg.  73]). 

Theorem  2  —  sufficient  conditions  for  convergence  of  the  Jacobi  and  Gauss- Seidel  relaxation 

Let  A  be  either  a  diagonally  dominant  or  an  irreducible  and  weakly  diagonally  dominant 
matrix.  Then,  both  the  associated  Jacobi  and  Gauss-Seidel  relaxation  methods  of  (4)  and  (5)  are 
convergent. 

Next,  we  turn  our  attention  to  the  successive  overrelaxation  method.  Since  the  inverse  of  a 
triangular  matrix  is  also  a  triangular  matrix,  and  its  determinant  is  equal  to  the  product  of  its 
diagonal  elements,  we  have 

det(Gw)  =  det[(I  —  wD-'LJ-'JdetKl  -  w)I  +  wD~lU]  =  (1  —  w)N. 

Since  det(G^)  =  fjfLi  it  follows  that 

N 

max  |X,|  >  |1  —  w|, 

t»i 

which,  by  theorem  1,  leads  to  the  following. 

Ihrorr-.  J  convergence  of  the  SOR  method 

p(G,J  >  |w-l| 


•  mt  mat  unreduced  bjr  Frobenius  for  matrices  which  (informally  speaking) 

.  .uni  wbosr  solut.ons  cannot  be  reduced  to  the  section  of  two  systems  of 
”^ducibi»  r^s'rices  when  Hisrrelising  boundary  value  problems  over 
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therefore  the  successive  overrelaxation  method  (6)  can  converge  only  if  0  <  u>  <  2. 

A  set  of  necessary  and  sufficient  conditions  for  convergence  for  the  successive  overrelaxation 
method  are  stated  in  the  following  theorem. 

Theorem  4  —  convergence  of  SOR  method  for  symmetric,  positive  definite  A 

If  A  is  real,  symmetric  matrix  with  positive  diagonal  elements,  then  the  successive  overrelaxation 
method  ( 5 )  is  convergent  if  and  only  if  0  <  w  <  2  and  A  is  positive  definite. 

The  same  conditions  for  convergence  hold  for  the  Gauss-Seidel  method  since,  by  definition, 
it  is  a  special  case  of  successive  overrelaxation. 

Corollary  2  —  convergence  of  Causs-Scidel  relaxation  for  symmetric,  positive  definite  A 

If  A  is  a  symmetric  matrix  with  positive  diagonal  elements,  then  the  Gai  ss-Seidel  method  is 
convergent  if  and  only  if  A  is  positive  definite. 

It  is  important  to  note  that  the  same  statement  cannot  be  made  about  the  Jacobi  method. 

B.2.  Basic  Gradient  Methods 

In  this  section  we  investigate  a  class  of  iterative  methods  which  are  naturally  associated 
with  optimization  theory.  These  are  the  so  called  gradient  methods  which,  in  their  full  generality, 
are  iterative  techniques  for  minimizing  nonlinear  functionals.  They  may  also  be  thought  of  as 
methods  for  solving  systems  of  linear  equations  for  the  special  case  where  the  functional  to  be 
minimized  is  a  quadratic  form. 

Assume  that  A  £  fflNN  is  symmetric.  Now  suppose  that  we  attempt  to  solve  the  following 
unconstrained  minimization  problem  involving  the  quadratic  form  £ (v) 

£(«)  =  £(v)  =  Av)  -  (f,  v), 

where  (•,  ■):  fStN  X  $lN  >-*  3?  denotes  the  familiar  Euclidean  inner  product.  From  the  well-known 
theorem  of  optimization  theory  (see,  e  g.  [Luenberger,  1973]),  the  gradient  vector  of  <f(u) 

V<f(u)  =  Au  -  f, 

vanishes  at  a  minimum  u  and,  moreover,  the  minimum  exists  and  is  unique  if  the  Hessian  matrix 

'£!£001  a. 

du,duj 

is  positive  definite. 

Thus,  for  symmetric,  positive  definite  A,  solving  the  minimization  problem  is  equivalent 
to  solving  the  system  of  linear  equations  Au  =  f  and,  consequently,  relaxation  methods  can  be 
thought  of  as  being  methods  for  descending  to  the  minimum  u  of  a  quadratic  functional. 
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Gradient  Descent 

Consider  the  iterative  method  =  u***  f  alfc*d*fc*,  where  at  each  iteration,  we  take 

a  step  in  the  direction  of  the  vector  d***.  To  minimize  £(u)  quickly,  we  should  move  in  the 
direction  of  steepest  descent,  which  is  given  by  the  negative  gradient;  that  is,  d***  =  —  V£'(u*k')  = 
f  —  Au***  =  r^k\  where  r*fc*  is  the  familiar  residual  vector.  Thus  we  obtain  the  nonatationary 
iterative  method  defined  by 


u(fc+0  =  u(*)  a(fc)r{fc). 


(10) 


It  seems  reasonable  to  choose  the  step  size  at  each  iteration  a*fc*  so  as  to  minimize  a(*)r(fc)). 

The  appropriate  value  can  easily  be  shown  to  be 


a 


(H*),  Ar<*>)‘ 


(11) 


On  the  other  hand,  if  we  fix  a**1*  =  a  for  all  iterations,  then  (10)  becomes 


u<*+‘>  =  (I  —  aA)U(t)  +  of,  (12) 

which  can  be  identified  as  a  stationary  iterative  method  with  iteration  matrix  G  =  I  —  oA. 

The  Conjugate  Gradient  Method 

Hestenes  and  Stiefel  [1952]  introduced  the  conjugate  gradient  method,  a  modification  of  the 
method  of  gradient  descent.  The  method  is  based  on  determining  vectors  d*0*,  d*1*, . . . ,  d*n—1* 
which  are  pairwise  conjugate  in  the  sense  that  (d(,\ Ad,J*)  =  0  for  i  ^  j.  The  ease  in  applying 
the  method  derives  from  the  fact  that  these  vectors  may  also  be  determined  iteratively.  With  the 
residual  vector  at  the  k ^  iteration  given  by  r***  —  t —  Au*k*  and  with  d*°*  =  r*0*,  the  algorithm 
for  determining  the  d(**  and  u***  is  as  follows. 


’.d1*1) 


d<* 


'Ad'*") 

d«*»  -  r<*>  - 


-i) 


o  <  k  <  N  -  1; 

1  <  *  <  N  -  1. 


Conj’ig.tcy  of  the  vectors  can  be  verified  along  with  the  fact  that  (r*1*,  r* 1  *)  =  0  for  t  ^  j.  This 
implies  that  r***  =  0  for  some  k  <  N .  Therefore,  in  the  absence  of  roundoff  errors,  the  method 
converges  in  at  most  N  iterations.  Of  course,  this  property  is  of  no  real  value  to  us  because  we 
must  deal  with  cases  where  N  is  very  large.  Nevertheless,  for  N  large,  typically  U1**  =s  u  for 
i<JV  and  the  algorithm  may  be  used  in  the  iterative  spirit.  The  conjugate  gradient  method  is 
certainly  the  most  expensive  of  the  algorithms  discussed  both  in  terms  of  space  (since  the  vectors 
„(*),  d *'  ),  r**',  and  Ad*1'*  must  be  stored)  and  in  terms  of  the  number  of  operations  to  complete 
one  iteration.  While  it  is  true  that,  for  model  problems,  the  number  of  iterations  required  to 
reduce  the  error  by  a  specified  amount  is  usually  considerably  less  than  for  the  other  methods,  the 
conjugate  gradient  method  seems  to  exceed  at  least  the  successive  overrelaxation  method  (with 
optimal  w)  in  total  number  of  operations  required  [see  e.g.  Young  and  Gregory,  1973,  pg.  1071). 
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Convergence  and  Comparisons  to  Relaxation 

It  is  easy  to  show  [Luenbcrger,  1973]  that  gradient  descent  (10)  with  (11)  is  convergent  for 
a  positive  definite  matrix  A  to  the  solution  u  =  A- ‘f  where  the  quadratic  form  £  is  minimised. 
On  the  other  hand,  we  must  be  a  little  more  careful  with  the  fixed-a  descent  algorithm  (12).  The 
eigenvalues  X,  of  its  iteration  matrix  G  are  related  to  the  eigenvalues  X]  of  A,  all  of  which  are 
positive  (since  A  is  positive  definite),  by  X,  =  1  —  aX[.  Therefore,  according  to  theorem  1,  we 
obtain  the  following. 

Theorem  5  —  necessary  and  sufficient  condition  for  convergence 

For  a  positive  definite  matrix  A,  the  fixed-a  method  of  equation  (IS)  is  convergent  if  and 
only  xfO  <  a  < 

Of  course,  convergence  is  quickest  for  that  a  which  minimizes  p(G ),  which  often  cannot  be 
determined  in  practice. 

By  comparing  (12)  with  the  Jacobi  method  (3),  we  can  convince  ourselves  that  the  two 
methods  become  identical  when  A  is  positive  definite,  has  identical  elements  on  its  main  diagonal 
(i.e.  D  =  al),  and  a  =  1/a.  Moreover,  Forsythe  and  Wasow  [1960,  pg.  239]  (see  also  |Milne, 
1970])  show  that  the  Gauss-Seidel  and  successive  overrelaxation  methods  are  also  subject  to 
interpretations  as  descent  methods.  In  these  cases,  however,  we  have  not  one,  but  a  set  of  direction 
vectors  which  turn  out  to  be  i„  the  coordinate  unit  vectors  of  SN .  During  each  iteration,  we 
take  a  sequence  of  steps  of  different  sizes  in  each  of  these  directions.  The  step  sizes  are  such  that 
f(u(‘)  -{-  qJ^'x,)  is  minimized,  and  are  given  by  oj**  —  r’*'/o„  for  the  Gauss-Seidel  method  and 
a]**  =  uir<l**/al,  for  the  successive  overrelaxation  method,  where  rjfc)  is  the  i1*1  element  of  the 
residual  vector  at  iteration  k. 

The  computation  of  an  optimal  a**)  at  each  iteration  according  to  (11)  requires  that 
be  stored  and  that  Ar,lt*  be  evaluated.  This  doubles  the  amount  of  work  required  per  iteration 
in  comparison  with  a  fixed-a  algorithm  or  Jacobi  relaxation.  This  raises  the  question:  which  is 
better  in  the  long  run;  N  iterations  of  gradient  descent,  or  2 N  iterations  of  the  fixed-a  or  Jacobi 
relaxation?  To  quantitatively  decide  the  issue,  a  convergence  analysis  ought  to  be  attempted.  This 
is  problem-dependent  and  is  generally  difficult  to  do,  so  we  will  simply  note  that  Forsythe  and 
Wasow  [1960,  pg.  225]  do  not  recommend  the  optimisation  of  a  and,  moreover,  refer  to  a  result  by 
Stiefel  [1955]  indicating  that  it  is,  at  best,  a  short-sighted  strategy.  Interestingly  enough,  Grimson 
[1981a]  used  the  optimal-a  algorithm  in  bis  implementation  of  the  surface  interpolation  algorithm 
(with  a  minor  modification  to  make  certain  that  the  fixed  constraints  are  never  modified,  thus, 
in  effect,  treating  them  as  essential  boundary  conditions).  Considering  the  statements  of  Forsythe 
and  Wasow,  it  is  not  surprising  that  extremely  slow  convergence  was  observed  in  spite  of  the 
extra  work  expended  at  each  iteration  to  compute  the  optimal  a. 
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C.  LOCAL  FOURIER  ANALYSIS  OF  RELAXATION 


In  this  appendix,  we  present  the  details  of  a  local  Fourier  analysis  of  the  Gauss- Seidel 
relaxation  and  obtain  the  smoothing  factor  for  this  method.  The  analysis  involves  studying 
separately  the  convergence  of  the  high-frequency  Fourier  components.  Since  these  components 
have  short  coupling  ranges,  we  can  perform  the  analysis  in  the  interior  of  ft,  ignoring  the  effects 
of  the  boundary  and  the  constraints. 

According  to  equation  (4.9),  the  minimising  displacement  u,^  at  an  interior  node  (», j)  is 
related  to  the  other  nodes  by 


20ut|J  —  8(ui_i,j  -f-  Ui+i.j  -(-  u,(j — i  +  u^j-j-j) 

+  2(u,_i0_i  -f  U|-f_] tj — i  4-  u,_i  lJ+i  -|-  u.+i.j  f  i) 
+  l(u«— 2,j  4“  4“  U,,j  —  2  4"  Ui,j+2)  =  0 


(1) 


(for  convenience,  we  have  suppressed  the  superscripts  h).  According  to  the  discussion  in  Appendix 
B,  at  iteration  k  of  the  Gausa-Seidel  relaxation  method,  v|^  is  replaced  by  a  new  value 
such  that 


20v«‘+»  -  8(vJirj  4-  +  tJ£.V  + 

4- 2(v['Lt!J)-i  4-  4-  4-v|^,,;+1) 

4-  +  v«+a.;  +  2  +  v!*H-a)  =  0 

The  errors  of  the  approximation  at  iteration  k  and  iteration  it  4-  1  are  given  by 


(2) 


.<*>  = 


u., 


and 


e,k+,)-u  _T,fc+,) 

e>,J 


respectively.  Subtracting  (2)  form  (1),  we  obtain 


9nJ'c+l)  j_J*4-i)  ,  „(*)  ^ 

20e..;  ~  8\ei-l.r  +ei-M.;  +  ei.;-l  4*  l ) 

4-  4-  4-  e,_,  fl  4-  e,  +  l  ;  + 

4- ,  .!*)  .  JM-i)  ,  <*)  \_0 

4^  ‘\et— a,j  4- +  ei,j— 2  o  ei,r+2 J  ~  °> 


(3) 


Suppose  that  the  error  consists  of  only  a  single  spatial  Fourier  component  <3  =  [u»j,  ws].  Then  the 
errors  at  node  (s',  j )  before  and  after  the  iteration  are  given  by 


respectively,  where  i  —  \J —■  1. 


(4) 


Substituting  (4)  into  (3),  dividing  through  by  e*(wi*+«».»)|  collecting  terms  pertaining  to 
the  same  iteration,  we  obtain 


84 


TERZOPOULOS 


MULTI-LEVEL  SURFACE  RECONSTRUCTION 


8(e"J'  -f  e1"’)  -f  2(e‘(u,,+u'j)  +  e*(— "•+"»>)  -(-  (e2‘"i  +  e2lu,)j 

+  v4a +,^20  —  -f  e~lWl)  +  2(e‘("1-w,)  -f  «*(-“!-“>))  -|-  (e— 2*"»  -f  e~ 2,WJ)^  =  0. 


The  amplification  of  the  u3  component  is  then  given  by 


Hgs(Q)  = 


4k+,) 

4k) 

8(e“*’1  -f-  e*-**’*)  —  2(e*l"1+Wj)  +  e‘l~ •*»+«* )J  —  (ei,Wl  +  e2*"*) 

20  —  8(e~ +  e— *“»)  +  2(e‘<u',— w*>  ■+■  e‘(—w>— "»))  +  (e— 2tw»  +  t~ 2tw>) 


Let  |b)|  =  max(|wj  |,  |w2|).  The  Gauss-Seidel  smoothing  factor  is  defined  as  the  smallest 
amplification  attained  for  a  high-frequency  component  of  the  error;  that  is,  a  component  which 
is  visible  on  the  fine  grid,  but  is  aliased  on  the  coarse  grid: 

Pas  =  „  JWi*. 

Evaluating  this  expression  numerically,  we  obtain  JiGS  ss  0.8  (for  u»i  =  1.6  and  w 2  =  0.3). 
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